using System;
using System.Collections.Generic;
using System.Windows.Forms;
using System.Drawing;
using System.Drawing.Drawing2D;
using System.IO;
using System.Runtime.Serialization;
using System.Runtime.Serialization.Formatters.Binary;
using System.Collections;
using System.Reflection;
using System.Security.Permissions;
using System.ComponentModel;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;
using System.Data;
using System.Data.Odbc;
using NumericalRecipes;
using FftLib;
using BmgNet;
using StaticUtility;

namespace ObjectLib
{
	[Serializable]
	public class Ribbon
	{
		public const uint MAX_PIX = 4095u;
		/// <summary>
		/// Holds a static resource for doing calculations in color space.
		/// </summary>
		private static KodakColorSpace staticKodakResource = new KodakColorSpace();
		public static ArrayList sectionFinder=new ArrayList();
		private static float[] rgb = new float[3];
		private static float[] logLms = new float[3];
		private static float[] lumByRg = new float[3];
		public string CategoryName;
		public string ProjName;
		public string PoName;
        public double PoDiameter;
		/// <summary>
		/// The points of the spine (center of the ribbon).
		/// </summary>
		public PointF[] Spine;
		/// <summary>
		/// The pixel locations.
		/// </summary>
		public PointF[,] Pixels;
		/// <summary>
		/// The pixel values of the ribbon, rgba values
		/// </summary>
		public float[,,] PixelValues;
		/// <summary>
		/// The number of pixels in the ribbon.
		/// </summary>
		public int PixelCt;
		/// <summary>
		/// The pixels that were included in the ribbon's rectangle.
		/// </summary>
		public bool[,] PixelIncluded;
		/// <summary>
		/// Gives the vector of directional distances that the ribbon will sample.  Units are pixels.  Positive values are
		/// inside of the leaf and negative values are outside of the leaf.  When a directional distance happens to violate
		/// this inside/outside relation, the pixel value is set to a null or empty value.
		/// The number of elements here determines the height (number of pixels) for the display ribbon.
		/// </summary>
		public float[] DirectionalDistances;
        /// <summary>
        /// The pixels with zero and negative directional distances.
        /// </summary>
        public PixelTagger PixNeg;
        /// <summary>
        /// The pixels with zero and positive directional distances.
        /// </summary>
        public PixelTagger PixPos;
        public double[] PixNeg_Mean,PixPos_Mean;
        public double[,] PixNeg_Cov,PixPos_Cov;
		/// <summary>
		/// Kernel for first-derivative calculation.  Gaussian over the indexes of Spine.
		/// </summary>
		public double SpineKernelSd;
		/// <summary>
		/// The space (in pixels) between successive points on the spine.  This should be the same number
		/// for all successive points.
		/// </summary>
		public double SpineSpacing;
		/// <summary>
		/// Holds the derivative of the edge (in radians).
		/// </summary>
		public double[] FirstAngularDerivative;
		/// <summary>
		/// Holds the curvature of the edge (in radians).
		/// </summary>
		public double[] SecondAngularDerivative;
		/// <summary>
		/// Provides Lum, Yb, and Rg values in the LMS-spaced derived from eigen-decomposition of natural images.
		/// </summary>
		public float[,] Lum, Yb, Rg;
		/// <summary>
		/// Gives the log-cone space like Ruderman, Cronin, and Ciao.
		/// </summary>
		public float[,] LogLum, LogYb, LogRg, Hue, Chroma, Saturation;
		/// <summary>
		/// The mean and standard deviation of the LMS values in color ribbons.
		/// </summary>
		public float MeanLum, SdLum, MeanRg, SdRg, MeanYb, SdYb;
		/// <summary>
		/// The mean and standard deviation of the log LMS values in color ribbons.
		/// </summary>
		public float MeanLogLum, SdLogLum, MeanLogRg, SdLogRg, MeanLogYb, SdLogYb;
		/// <summary>
		/// The mean and standard deviation of the log LMS values in color ribbons.
		/// </summary>
		public float MeanChroma, SdChroma, MeanSat, SdSat;
		/// <summary>
		/// The categorization values for this Ribbon.
		/// </summary>
		public object CatData;
        /// <summary>
        /// The filter responses for each (1) filter type, (2) color channel, (3) position along the ribbon.
        /// </summary>
        public ComplexD[][][] Filters;
        /// <summary>
        /// Evaluates a filter at the specified coordinate (y, x) from center as if the filter was first rotated at the specified angle.
        /// </summary>
        /// <param name="filter">A filter.</param>
        /// <param name="y">The offset from the center of the filter in the first index (vertical spatial direction)</param>
        /// <param name="x">The offset from the center of the filter in the second index (horizontal spatial direction)</param>
        /// <param name="angle">The rotation of the filter.</param>
        /// <returns>The filter value (nearest neighbor interpolation).</returns>
        public static ComplexD Evaluate(ComplexD[,] filter,double y,double x,double angle)
        {
            double c = Math.Cos(angle);
            double s = Math.Sin(angle);
            double xx = c * x - s * y; //  rotation is backwards, as the patch rotates one way, the point rotates the other.
            double yy = s * x + c * y;
            int ni=filter.GetLength(0),nj=filter.GetLength(1);

            double iy = yy + (double)(ni - 1) / 2.0;
            double jx = xx + (double)(nj - 1) / 2.0;
            if(iy < 0.0 || iy >= (double)(ni - 1) || jx < 0.0 || jx >= (double)(nj - 1))
                return ComplexD.Zero;

            int i = (int)iy;
            int j = (int)jx;

            double AB,aB,ab,Ab,dy,dx;
            dy = iy % 1.0;
            dx = jx % 1.0;
            AB = (1.0 - dy) * (1.0 - dx);
            Ab = (1.0 - dy) * dx;
            aB = dy * (1.0 - dx);
            ab = dy * dx;
            return AB * filter[i,j] + Ab * filter[i + 1,j] + aB * filter[i,j + 1] + ab * filter[i + 1,j + 1];
        }
        /// <summary>
        /// A cache for the window evaluation, used in the 'Evaluate' function.
        /// </summary>
        private static ComplexD[,] _winCache;
        /// <summary>
        /// Computes the convolution of a filter and an image at a specified point in LogLms-Opponent space.  Assumes the filter is super-sampled
        /// (at least twice as large width and height).
        /// </summary>
        /// <param name="window">The filter to convolve against the image.</param>
        /// <param name="rgbFac">The multiplication factor for RGB pixels.</param>
        /// <param name="pt">PointF  The image point for convolution.</param>
        /// <param name="len">int  The diameter of the filter (usually smaller than filter.Ni or filter.Nj, facilitates super-sampling).</param>
        /// <param name="angle">double The filter is rotated at this angle (in radians).</param>
        /// <returns>The convolution of a filter on an image for each color channel of the image.  The real values are returned in the first column, imaginary in second column.</returns>
        public static ComplexD[] Convolve(Ushort3 pix,double rgbFac,PointF pt,int len,double angle,ComplexD[,] filter)
        {
            int i,j,ii,jj,pi,pj;
            int ni = filter.GetLength(0);
            int nj = filter.GetLength(1);
            float[] temp = new float[3],dum = new float[3],lab = new float[3];
            ComplexD[] output = new ComplexD[3];
            ComplexD win;
            ComplexD sum = ComplexD.Zero,mean;
            ComplexD over = ComplexD.Zero,under = ComplexD.Zero,area;
            if(_winCache == null || _winCache.GetLength(0) != len)
                _winCache = new ComplexD[len,len];

            double y,x;

            //  shift from the center of the filter to upper left
            pt.X -= (float)(len-1) / 2.0f;
            pt.Y -= (float)(len-1) / 2.0f;

            pj = (int)pt.X;  //  The point on the image pj + a, is where the patch's left side is.
            pi = (int)pt.Y;  //  The point on the image pi + b, is where the patch's top side is.
            float a = pt.X % 1.0f;
            float b = pt.Y % 1.0f;

            float ab = a * b;
            float Ab = (1.0f - a) * b;
            float aB = a * (1.0f - b);
            float AB = (1.0f - a) * (1.0f - b);

            int ct = 0;

            for(i = 0;i < len;i++)
            {
                //  Get the unrotated filter coordinate.
                y = (double)i / (double)(len - 1) - 0.5;
                for(j = 0;j < len;j++)
                {
                    //  Get the unrotated filter coordinate.
                    x = (double)j / (double)(len - 1) - 0.5;
                    if(y * y + x * x > 0.5)
                        _winCache[i,j] = ComplexD.Empty;
                    else
                    {
                        ct++;
                        //  Get the filter value
                        win = Evaluate(filter,y * (double)ni,x * (double)nj,angle);
                        sum += win;
                        if(win.Real > 0.0)
                            over.Real += win.Real;
                        else
                            under.Real += win.Real;
                        if(win.Imag > 0.0)
                            over.Imag += win.Imag;
                        else
                            under.Imag += win.Imag;
                        _winCache[i,j] = win;
                    }
                }
            }
            mean = sum / (double)(ct*ct);
            area = over - under;
            for(i = 0;i < len;i++)
            {
                for(j = 0;j < len;j++)
                {
                    if(_winCache[i,j] != ComplexD.Empty)
                    {
                        //  For non-zero filter values only
                        if(_winCache[i,j] != ComplexD.Zero)
                        {
                            //	(0,0)
                            ii = pi + i;    //  ii-->A, ii+1-->a
                            jj = pj + j;    //  jj-->B, jj+1-->b

                            if(ii >= 0 && jj >= 0 && ii + 1 < pix.Ni && jj + 1 < pix.Nj
                                && pix[ii,jj,0] < 4095u && pix[ii,jj,1] < 4095u && pix[ii,jj,2] < 4095u
                                && pix[ii + 1,jj,0] < 4095u && pix[ii + 1,jj,1] < 4095u && pix[ii + 1,jj,2] < 4095u
                                && pix[ii,jj + 1,0] < 4095u && pix[ii,jj + 1,1] < 4095u && pix[ii,jj + 1,2] < 4095u
                                && pix[ii + 1,jj + 1,0] < 4095u && pix[ii + 1,jj + 1,1] < 4095u && pix[ii + 1,jj + 1,2] < 4095u)
                            {
                                //  Adjust filter
                                win = (_winCache[i,j] - mean);
                                win.Real /= area.Real;
                                win.Imag /= area.Imag;

                                //  Interpolate
                                //  (0,0)
                                dum[0] = (float)(rgbFac * (double)pix[ii,jj,0]);
                                dum[1] = (float)(rgbFac * (double)pix[ii,jj,1]);
                                dum[2] = (float)(rgbFac * (double)pix[ii,jj,2]);
                                KodakColorSpace.RgbToLogLms(dum,temp);
                                KodakColorSpace.LogLmsToLumYbRg(temp,dum);

                                lab[0] = AB * dum[0];
                                lab[1] = AB * dum[1];
                                lab[2] = AB * dum[2];

                                //	(1,0)
                                ii++;
                                dum[0] = (float)(rgbFac * (double)pix[ii,jj,0]);
                                dum[1] = (float)(rgbFac * (double)pix[ii,jj,1]);
                                dum[2] = (float)(rgbFac * (double)pix[ii,jj,2]);
                                KodakColorSpace.RgbToLogLms(dum,temp);
                                KodakColorSpace.LogLmsToLumYbRg(temp,dum);

                                lab[0] += aB * dum[0];
                                lab[1] += aB * dum[1];
                                lab[2] += aB * dum[2];

                                //	(1,1)
                                jj++;
                                dum[0] = (float)(rgbFac * (double)pix[ii,jj,0]);
                                dum[1] = (float)(rgbFac * (double)pix[ii,jj,1]);
                                dum[2] = (float)(rgbFac * (double)pix[ii,jj,2]);
                                KodakColorSpace.RgbToLogLms(dum,temp);
                                KodakColorSpace.LogLmsToLumYbRg(temp,dum);
                                lab[0] += ab * dum[0];
                                lab[1] += ab * dum[1];
                                lab[2] += ab * dum[2];

                                //	(0,1)
                                ii--;
                                dum[0] = (float)(rgbFac * (double)pix[ii,jj,0]);
                                dum[1] = (float)(rgbFac * (double)pix[ii,jj,1]);
                                dum[2] = (float)(rgbFac * (double)pix[ii,jj,2]);
                                KodakColorSpace.RgbToLogLms(dum,temp);
                                KodakColorSpace.LogLmsToLumYbRg(temp,dum);
                                lab[0] += Ab * dum[0];
                                lab[1] += Ab * dum[1];
                                lab[2] += Ab * dum[2];

                                //	Sum
                                output[0] += (double)lab[0] * win;
                                output[1] += (double)lab[1] * win;
                                output[2] += (double)lab[2] * win;
                            }
                        }
                    }
                }
            }
            return output;
        }
		/// <summary>
		/// Constructs a ribbon.
		/// </summary>
		/// <param name="po">The parsed object associated with this ribbon.</param>
		/// <param name="image">float[,,]  Image pixels.</param>
		/// <param name="Spine">PointF[]  The polygon spine.</param>
		/// <param name="DirectionalDistances">float[]  The directional distance values of the ribbon.  Negative
		/// outside polygon, positive inside.  Should go from negative to positive, monotonically, with a 'zero' in the middle.</param>
		/// <param name="SpineKernelSd">double  The standard deviation of the smoothing kernel for the spine's derivative values.</param>
        /// <param name="filter">Complex2f  The filter (usually log-gabor) with even (real) and odd (imaginary) phases.</param>
        /// <param name="filterLength">int  The length of the filter (usually smaller than filter.Ni or filter.Nj, facilitates super-sampling).</param>
        /// <param name="filterOffsets">PointF  Each filter is offset by this amount in pixels, where filterOffsets[.].X is offset
        /// parallel to the contour and filterOffsets[.].Y is offset perpendicular to the contour</param>
		public Ribbon( ParsedObject po, PathF pf, Ushort3 corePixels, float rgbFac, double SpineSpacing,
			PointF[] Spine, float[] DirectionalDistances, double SpineKernelSd,
            ComplexD[][,] filters, int filterLength, PointF[] filterOffsets )
		{
            int i,j,ix,iy,index,k,jump;

			//	The Spine and DirectionalDistances define this ribbon.
			this.Spine = Spine;
			this.DirectionalDistances = DirectionalDistances;
			this.SpineSpacing = SpineSpacing;

			//	The object name
			PoName = po.Name;
            PoDiameter = Math.Sqrt((double)po.Area);
			ProjName = po.Category.Owner.ShortProjname;
			CategoryName = po.Category.Name;

			//	The categorization values
			if( pf.GetType()==typeof(ContourF) )
				CatData = ((ContourF)pf).CatData;
			else
				CatData = po.CatData;

			//	Allocate memory for the interpolated pixels
			Pixels = new PointF[DirectionalDistances.Length, Spine.Length];
			PixelValues = new float[DirectionalDistances.Length, Spine.Length, 4];
			PixelIncluded = new bool[DirectionalDistances.Length, Spine.Length];

			//	Allocate memory for the derivatives
			FirstAngularDerivative = new double[Spine.Length];
			SecondAngularDerivative = new double[Spine.Length];

            //  Allocate memroy for the pixel taggers
            PixNeg = new PixelTagger();
            PixPos = new PixelTagger();

            int iLim = corePixels.Ni-1;
			int jLim = corePixels.Nj-1;
			int ribbonHeight = DirectionalDistances.Length;
			PixelCt=0;
			double theta;
			float cos, sin, dx, dy, dxx, dyy;
			bool inPoly=false;
			PointF pt=new PointF(0.0f,0.0f);
			//	for each point along the Spine
			for(i=0; i<Spine.Length; i++)
			{
				//	get the angle of the Spine, bound within [-pi, pi]
				theta = Util.KernelDerivative(Spine, SpineKernelSd, i, false);
				if(theta > Math.PI)
					theta-=2.0*Math.PI;
				else if(theta<-Math.PI)
					theta+=2.0*Math.PI;
				FirstAngularDerivative[i] = theta;
				//	define the angle of the perpendicular
				theta-=Math.PI/2.0;
				cos=(float)Math.Cos(theta);
				sin=(float)Math.Sin(theta);
				//	for each sample along the perpendicular
				for(j=0; j<ribbonHeight; j++)
				{
					Pixels[ribbonHeight-1-j,i]
						=pt
						=new PointF(
						Spine[i].X+DirectionalDistances[j]*cos,
						Spine[i].Y+DirectionalDistances[j]*sin );
					//	If pixels are only included by object, see if this object contains the pixel.
					if( pf.GetType()==typeof(PolygonF) )
						inPoly = ((PolygonF)pf).InPoly( pt.X, pt.Y ) > 0;
					else
						inPoly = ((ContourF)pf).PolyOwner.InPoly( pt.X, pt.Y ) > 0;
					//	insure positive DirectionalDistances are positive inside the polygon and negative outside.
					if(  (pf.GetType()==typeof(ContourF) && !inPoly)
						|| pf.GetType()==typeof(PolygonF) && 
							(  (DirectionalDistances[j]>0.0f && inPoly)
							|| (DirectionalDistances[j]<0.0f && !inPoly) ) )
					{
						PixelValues[ribbonHeight-1-j,i,0]
							=PixelValues[ribbonHeight-1-j,i,1]
							=PixelValues[ribbonHeight-1-j,i,2]
							=PixelValues[ribbonHeight-1-j,i,3]=0.0f;
						PixelIncluded[ribbonHeight-1-j,i]=false;
					}
					else 
					{
						PixelIncluded[ribbonHeight-1-j,i]=true;
						//	2d interpolation
						ix = (int)pt.X;
						iy = (int)pt.Y;
						dx = pt.X-(float)ix;
						dy = pt.Y-(float)iy;
						dxx=1.0f-dx;
						dyy=1.0f-dy;
						index = ribbonHeight-1-j;

						try
						{
							if( iy >= iLim || ix >= jLim || iy<0 || ix<0
								||	corePixels[iy,ix,0]>=MAX_PIX
								||	corePixels[iy,ix,1]>=MAX_PIX
								||	corePixels[iy,ix,2]>=MAX_PIX
								||	corePixels[iy+1,ix,0]>=MAX_PIX
								||	corePixels[iy+1,ix,1]>=MAX_PIX
								||	corePixels[iy+1,ix,2]>=MAX_PIX
								||	corePixels[iy+1,ix+1,0]>=MAX_PIX
								||	corePixels[iy+1,ix+1,1]>=MAX_PIX
								||	corePixels[iy+1,ix+1,2]>=MAX_PIX
								||	corePixels[iy,ix+1,0]>=MAX_PIX
								||	corePixels[iy,ix+1,1]>=MAX_PIX
								||	corePixels[iy,ix+1,2]>=MAX_PIX )
							{
								PixelValues[index,i,0]
									=PixelValues[index,i,1]
									=PixelValues[index,i,2]
									=PixelValues[index,i,3]=0.0f;
								PixelIncluded[index,i]=false;
							}
							else
							{
								PixelValues[index,i,0]=
									dxx*dyy*(float)corePixels[iy,ix,0]*rgbFac
									+dxx*dy*(float)corePixels[iy,ix+1,0]*rgbFac
									+dx*dyy*(float)corePixels[iy+1,ix,0]*rgbFac
									+dx*dy*(float)corePixels[iy+1,ix+1,0]*rgbFac;
								PixelValues[index,i,1]=
									dxx*dyy*(float)corePixels[iy,ix,1]*rgbFac
									+dxx*dy*(float)corePixels[iy,ix+1,1]*rgbFac
									+dx*dyy*(float)corePixels[iy+1,ix,1]*rgbFac
									+dx*dy*(float)corePixels[iy+1,ix+1,1]*rgbFac;
								PixelValues[index,i,2]=
									dxx*dyy*(float)corePixels[iy,ix,2]*rgbFac
									+dxx*dy*(float)corePixels[iy,ix+1,2]*rgbFac
									+dx*dyy*(float)corePixels[iy+1,ix,2]*rgbFac
									+dx*dy*(float)corePixels[iy+1,ix+1,2]*rgbFac;
								PixelValues[index,i,3]=
									dxx*dyy*(float)corePixels[iy,ix,3]*rgbFac
									+dxx*dy*(float)corePixels[iy,ix+1,3]*rgbFac
									+dx*dyy*(float)corePixels[iy+1,ix,3]*rgbFac
									+dx*dy*(float)corePixels[iy+1,ix+1,3]*rgbFac;

								PixelCt++;
							}
						}
						catch
						{
							PixelValues[index,i,0]
								=PixelValues[index,i,1]
								=PixelValues[index,i,2]
								=PixelValues[index,i,3]=0.0f;
							PixelIncluded[index,i]=false;
						}
                        //  Add to pixel tagger
                        if( PixelIncluded[index,i] )
                        {
                            if( DirectionalDistances[index] < 0.0f )
                                PixNeg.AddPixel( (int)(Pixels[index,i].X+0.5f), (int)(Pixels[index,i].Y+0.5f) );
                            else if(DirectionalDistances[index] > 0.0f)
                                PixPos.AddPixel((int)(Pixels[index,i].X+0.5f),(int)(Pixels[index,i].Y+0.5f));
                            //  Add in-between pixels
                            if( i>0 && PixelIncluded[index,i-1] )
                            {
                                dx = Pixels[index,i].X - Pixels[index,i-1].X;
                                dy = Pixels[index,i].Y - Pixels[index,i-1].Y;
                                dxx = dx*dx+dy*dy;
                                if( dxx > 1.0 )
                                {
                                    dxx = (float)Math.Sqrt(dxx);
                                    dx/=dxx;
                                    dy/=dxx;
                                    if(DirectionalDistances[index] < 0.0f)
                                        for(k=1; k<(int)dxx; k++)
                                            PixNeg.AddPixel((int)(Pixels[index,i].X-dx*(float)k+0.5f),(int)(Pixels[index,i].Y-dy*(float)k+0.5f));
                                    else if(DirectionalDistances[index] > 0.0f)
                                        for(k=1;k<(int)dxx;k++)
                                            PixPos.AddPixel((int)(Pixels[index,i].X-dx*(float)k+0.5f),(int)(Pixels[index,i].Y-dy*(float)k+0.5f));
                                }
                            }
                        }
					}
				}
			}

			//	Now, guarantee that each point on the ribbon is no less than the directional distance from other points on the spine.
            //  Also, compute the surrounding polygons
			if(Spine.Length>2)
			{
				//	'k'	refers to a level of DirectionalDistances (vertical)			
				//	'i' refers to a spot along the tape (horizontal)
				//	'j' refers to another spot along the spine (horizontal)
				//	The loop enumerates all pairs (i,j) where i!=j
				PointF pfi, pfj;
				double sqdist,mindist,dist;
				for(k=0; k<DirectionalDistances.Length; k++)
				{
					mindist = DirectionalDistances[k];
					for(i=0; i<Spine.Length; i++)
					{
						if(PixelIncluded[k,i])
						{
							pfi = Pixels[k,i];
							for(j=0; j<Spine.Length; j++)
							{
								if(i!=j)
								{
									pfj = Spine[j];
									dx = pfj.X-pfi.X;
									dy = pfj.Y-pfi.Y;
									sqdist = dx*dx+dy*dy;
									dist = Math.Sqrt(sqdist);
									//	Pixels[k,i] is closer to Spine[j] than Spine[i], remove it
									//	from consideration.
									if(dist<mindist)
										PixelIncluded[k,i] = false;
									else
									{
										//	All susubsequent spine points within (dist-mindist)/SpineSpacing cannot be closer either.
										jump = (int)( (dist-mindist)/SpineSpacing )-1;
										if(jump>0)
											j+=jump;
									}
								}
							}
						}
					}
                }

                #region Positive and Negative PixelTagger
                //  Now compute mean and cov of each PixelTagger, checking for bad pixels and removing them from each PixelTagger.
                jump=0;
                PixelTagger pix=null;
                double[] mean;
                double[,] cov;
                while(jump<2)
                {
                    if(jump==0)
                        pix=PixPos;
                    else if(jump==1)
                        pix=PixNeg;
                    if(pix!=null && pix.NumPts > 10)
                    {
                        //  Initialize result
                        mean = new double[3];
                        cov = new double[3,3];
                        for(i=pix.Left;i<=pix.Right;i++)
                            for(j=pix.Top;j<=pix.Bottom;j++)
                                if(pix.PixelIsTagged(i,j))
                                {
                                    if(corePixels[j,i,0]>=4095 || corePixels[j,i,0]>=4095 || corePixels[j,i,0]>=4095)
                                        pix.RemovePixel(i,j);
                                    else
                                    {
                                        rgb[0] = corePixels[j,i,0]*rgbFac;
                                        rgb[1] = corePixels[j,i,1]*rgbFac;
                                        rgb[2] = corePixels[j,i,2]*rgbFac;
                                        KodakColorSpace.RgbToLogLms(rgb,logLms);
                                        KodakColorSpace.LogLmsToLumYbRg(logLms,lumByRg);
                                        mean[0]+=(double)lumByRg[0];
                                        mean[1]+=(double)lumByRg[1];
                                        mean[2]+=(double)lumByRg[2];
                                    }
                                }
                        if(pix.NumPts > 10)
                        {
                            mean[0] = mean[0]/(double)pix.NumPts;
                            mean[1] = mean[1]/(double)pix.NumPts;
                            mean[2] = mean[2]/(double)pix.NumPts;

                            for(i=pix.Left;i<=pix.Right;i++)
                                for(j=pix.Top;j<=pix.Bottom;j++)
                                    if(pix.PixelIsTagged(i,j))
                                    {
                                        rgb[0] = corePixels[j,i,0]*rgbFac;
                                        rgb[1] = corePixels[j,i,1]*rgbFac;
                                        rgb[2] = corePixels[j,i,2]*rgbFac;
                                        KodakColorSpace.RgbToLogLms(rgb,logLms);
                                        KodakColorSpace.LogLmsToLumYbRg(logLms,lumByRg);
                                        rgb[0]=lumByRg[0]-(float)mean[0];
                                        rgb[1]=lumByRg[1]-(float)mean[1];
                                        rgb[2]=lumByRg[2]-(float)mean[2];
                                        cov[0,0] += rgb[0]*rgb[0];
                                        cov[1,1] += rgb[1]*rgb[1];
                                        cov[2,2] += rgb[2]*rgb[2];
                                        cov[0,1] += rgb[0]*rgb[1];
                                        cov[0,2] += rgb[0]*rgb[2];
                                        cov[1,2] += rgb[1]*rgb[2];
                                    }

                            cov[0,0] /= (double)(pix.NumPts-1);
                            cov[1,1] /= (double)(pix.NumPts-1);
                            cov[2,2] /= (double)(pix.NumPts-1);
                            cov[0,1] /= (double)(pix.NumPts-1);
                            cov[0,2] /= (double)(pix.NumPts-1);
                            cov[1,2] /= (double)(pix.NumPts-1);
                            cov[1,0] = cov[0,1];
                            cov[2,0] = cov[0,2];
                            cov[2,1] = cov[1,2];

                            //  Store result
                            if(jump==0)
                            {
                                PixPos_Mean=mean;
                                PixPos_Cov=cov;
                            }
                            else if(jump==1)
                            {
                                PixNeg_Mean=mean;
                                PixNeg_Cov=cov;
                            }
                            else
                                throw new ApplicationException("error");
                        }
                        else
                            pix = null;
                    }
                    jump++;
                }
                #endregion
            }

			//	Compute the second angular derivative, or curvature; insure between -pi and pi.
			Util.AngularDiff(FirstAngularDerivative, SecondAngularDerivative, FirstAngularDerivative.Length, false);
            //  Compute filters.
            
            if(filters != null)
            {
                Filters = new ComplexD[filters.Length][][];
                for(i=0; i<filters.Length; i++)
                {
                    Filters[i] = new ComplexD[][] {
                        new ComplexD[Spine.Length],
                        new ComplexD[Spine.Length],
                        new ComplexD[Spine.Length] };
                }
                ComplexD[] filt;
                double angle;
                PointF spinePt, dev=new PointF(0.0f,0.0f);
                for(i=0; i<this.Spine.Length; i++)
                {
                    for(j=0; j<filters.Length; j++)
                    {
                        spinePt = this.Spine[i];
                        //  parallel angle
                        angle = FirstAngularDerivative[i];
                        dev.X = (float)Math.Cos(angle)*filterOffsets[j].X;
                        dev.Y = (float)Math.Sin(angle)*filterOffsets[j].X;
                        //  perpendicular angle
                        angle+=Util.HalfPi;
                        dev.X += (float)Math.Cos(angle)*filterOffsets[j].Y;
                        dev.Y += (float)Math.Sin(angle)*filterOffsets[j].Y;
                        //  offset
                        spinePt.X += dev.X;
                        spinePt.Y += dev.Y;
                        //  convolve
                        //  The sign of the first angular derivative is reversed because the image is upside-down compared to the filter.
                        filt = Convolve(corePixels, rgbFac, spinePt, filterLength, -FirstAngularDerivative[i]+Util.HalfPi, filters[j]);
                        Filters[j][0][i] = filt[0];
                        Filters[j][1][i] = filt[1];
                        Filters[j][2][i] = filt[2];
                    }
                }
            }
		}
		/// <summary>
		/// For Space==ColorSpace.LmsOpp
		/// (1) Calculates the luminance ribbon, Lum, LumMean, LumSd, and Painter.LumBmp
		/// (2) Calculates the yellow-blue ribbon, Yb, YbMean, YbSd, and Painter.YbBmp
		/// (3) Calculates the red-green ribbon, Rg, RgMean, RgSd, and Painter.RgBmp
		/// For Space==ColorSpace.LogLmsOpp
		/// (1) Calculates the log luminance ribbon, LogLum, LogLumMean, LogLumSd, and Painter.LogLumBmp
		/// (2) Calculates the log yellow-blue ribbon, LogYb, LogYbMean, LogYbSd, and Painter.LogYbBmp
		/// (3) Calculates the log red-green ribbon, LogRg, LogRgMean, LogRgSd, and Painter.LogRgBmp
		/// </summary>
		public void CalcColorRibbons(ColorSpace space)
		{
			//	Get ribbon dimensions.
			int rows=PixelValues.GetLength(0);
			int cols=PixelValues.GetLength(1);

			float N = (float)(PixelCt);

			if(space==ColorSpace.LmsOpp)
			{

				//	Allocate memory
				Lum = new float[rows,cols];
				Yb = new float[rows,cols];
				Rg = new float[rows,cols];

				//	Initialize means and standard deviations.
				MeanLum=SdLum=0.0f;
				MeanYb=SdYb=0.0f;
				MeanRg=SdRg=0.0f;

				int i,j;

				//	Fill ribbons and start computing the means
				for(i=0; i<rows; i++)
				{
					for(j=0; j<cols; j++)
					{
						if(this.PixelIncluded[i,j])
						{
							rgb[0] = PixelValues[i,j,0];
							rgb[1] = PixelValues[i,j,1];
							rgb[2] = PixelValues[i,j,2];

							KodakColorSpace.RgbToLms(rgb, logLms);
							KodakColorSpace.LmsToLumYbRg(logLms, lumByRg);

							Lum[i,j]=lumByRg[0];
							MeanLum+=lumByRg[0];

							Yb[i,j]=lumByRg[1];
							MeanYb+=lumByRg[1];

							Rg[i,j]=lumByRg[2];
							MeanRg+=lumByRg[2];
						}
					}
				}

				//	Finish computing means
				MeanLum/=N;
				MeanYb/=N;
				MeanRg/=N;

				double diff;
				double ssLum=0.0, ssYb=0.0, ssRg=0.0;

				//	Compute the sum-squared error
				for(i=0; i<rows; i++)
				{
					for(j=0; j<cols; j++)
					{
						if(this.PixelIncluded[i,j])
						{
							diff = Lum[i,j]-MeanLum;
							ssLum+=diff*diff;

							diff = Yb[i,j]-MeanYb;
							ssYb+=diff*diff;

							diff = Rg[i,j]-MeanRg;
							ssRg+=diff*diff;
						}
					}
				}

				//	Compute the standard deviations
				SdLum = (float)Math.Sqrt(ssLum/((double)N-1));
				SdYb = (float)Math.Sqrt(ssYb/((double)N-1));
				SdRg = (float)Math.Sqrt(ssRg/((double)N-1));
			}
			else if(space==ColorSpace.LogLmsOpp)
			{

				//	Allocate memory
				LogLum = new float[rows,cols];
				LogYb = new float[rows,cols];
				LogRg = new float[rows,cols];
				Hue = new float[rows,cols];
				Chroma = new float[rows,cols];
				Saturation = new float[rows,cols];

				//	Initialize means and standard deviations.
				MeanLogLum=SdLogLum=0.0f;
				MeanLogYb=SdLogYb=0.0f;
				MeanLogRg=SdLogRg=0.0f;

				int i,j;

				//	Fill ribbons and start computing the means
				for(i=0; i<rows; i++)
				{
					for(j=0; j<cols; j++)
					{
						if(this.PixelIncluded[i,j])
						{
							rgb[0] = PixelValues[i,j,0];
							rgb[1] = PixelValues[i,j,1];
							rgb[2] = PixelValues[i,j,2];

							KodakColorSpace.RgbToLogLms(rgb, logLms);
							KodakColorSpace.LogLmsToLumYbRg(logLms, lumByRg);

							LogLum[i,j]=lumByRg[0];
							MeanLogLum+=lumByRg[0];

							LogYb[i,j]=lumByRg[1];
							MeanLogYb+=lumByRg[1];

							LogRg[i,j]=lumByRg[2];
							MeanLogRg+=lumByRg[2];

							Hue[i,j]=KodakColorSpace.LumYbRgToHue(lumByRg);

							Chroma[i,j]=KodakColorSpace.LumYbRgToChroma(lumByRg);
							MeanChroma+=Chroma[i,j];

							Saturation[i,j]=KodakColorSpace.LumYbRgToSat(lumByRg);
							MeanSat+=Saturation[i,j];
						}
					}
				}

				//	Finish computing means
				MeanLogLum/=N;
				MeanLogYb/=N;
				MeanLogRg/=N;
				MeanChroma/=N;
				MeanSat/=N;

				double diff;
				double ssLum=0.0, ssYb=0.0, ssRg=0.0, ssChroma=0.0, ssSat=0.0;

				//	Compute the sum-squared error
				for(i=0; i<rows; i++)
				{
					for(j=0; j<cols; j++)
					{
						if(this.PixelIncluded[i,j])
						{
							diff = LogLum[i,j]-MeanLogLum;
							ssLum+=diff*diff;

							diff = LogYb[i,j]-MeanLogYb;
							ssYb+=diff*diff;

							diff = LogRg[i,j]-MeanLogRg;
							ssRg+=diff*diff;

							diff = Chroma[i,j]-MeanChroma;
							ssChroma+=diff*diff;

							diff = Saturation[i,j]-MeanSat;
							ssSat+=diff*diff;
						}
					}
				}

				//	Compute the standard deviations
				SdLogLum = (float)Math.Sqrt(ssLum/((double)N-1));
				SdLogYb = (float)Math.Sqrt(ssYb/((double)N-1));
				SdLogRg = (float)Math.Sqrt(ssRg/((double)N-1));
				SdChroma = (float)Math.Sqrt(ssChroma/((double)N-1));
				SdSat = (float)Math.Sqrt(ssSat/((double)N-1));
			}
		}
		/// <summary>
		/// Demarcates the valid sections by each section's first and last position.  The section
		/// is valid if all pixels are included for all indexes into this.DirectionalDistances.
		/// A position is simiply an index into this.Spine.
		/// </summary>
		/// <param name="minimumLength">Specifies the minimum length of a valid section.</param>
		/// <returns>An array of valid section, where each section is a pair of integers with
		/// a given first and last position.  A position is simply an index into this.Spine.</returns>
		public int[][] ValidSections(int minimumLength)
		{
			sectionFinder.Clear();
			int validStart=-1;
			int width = (int)(this.RibbonWidth+0.5f);
			int i;
			for(i=0; i<width; i++)
			{
				if(ValidPosition(i))
				{
					if(validStart==-1)
						validStart=i;
				}
				else
				{
					//	Add the new section if it exceeds minimum length
					if(validStart>=0 && i-validStart>=minimumLength)
						sectionFinder.Add(new int[] {validStart, i-1});
					validStart=-1;
				}
			}
			//	Add the new section if it exceeds minimum length
			if(validStart>=0 && i-validStart>=minimumLength)
				sectionFinder.Add(new int[] {validStart, width-1});
			return (int[][])sectionFinder.ToArray(typeof(int[]));
		}
		/// <summary>
		/// Demarcates the valid sections by each section's first and last position.  The section
		/// is valid if all pixels are included for indexes into this.DirectionalDistances of
		/// all ind1 and ind2.  A position is simiply an index into this.Spine.
		/// </summary>
		/// <param name="ind1">Indexes into DirectionalDistances.</param>
		/// <param name="ind2">Indexes into DirectionalDistances.</param>
		/// <param name="minimumLength">Specifies the minimum length of a valid section.</param>
		/// <returns>An array of valid section, where each section is a pair of integers with
		/// a given first and last position.  A position is simply an index into this.Spine.</returns>
		public int[][] ValidSections(int[] ind1, int[] ind2, int minimumLength)
		{
			sectionFinder.Clear();
			int validStart=-1;
			int width = (int)(this.RibbonWidth+0.5f);
			int i;
			for(i=0; i<width; i++)
			{
				if(ValidPosition(i,ind1,ind2))
				{
					if(validStart==-1)
						validStart=i;
				}
				else
				{
					//	Add the new section if it exceeds minimum length
					if(validStart>=0 && i-validStart>=minimumLength)
						sectionFinder.Add(new int[] {validStart, i-1});
					validStart=-1;
				}
			}
			//	Add the new section if it exceeds minimum length
			if(validStart>=0 && i-validStart>=minimumLength)
				sectionFinder.Add(new int[] {validStart, width-1});
			return (int[][])sectionFinder.ToArray(typeof(int[]));
		}
		/// <summary>
		/// Determines whether the position (index into this.Spine) has valid entries
		/// for all indices into this.DirectionalDistances.
		/// </summary>
		/// <param name="position"></param>
		/// <returns></returns>
		public bool ValidPosition(int position)
		{
			int i;
			for(i=0; i<this.DirectionalDistances.Length; i++)
				if(!this.PixelIncluded[i,position])
					return false;
			return true;
		}
		/// <summary>
		/// Determines whether the position (index into this.Spine) has valid entries for all indices into
		/// this.DirectionalDistances enumerated in ind1 and ind2 together.
		/// </summary>
		/// <param name="position"></param>
		/// <param name="ind1"></param>
		/// <param name="ind2"></param>
		/// <returns></returns>
		public bool ValidPosition(int position, int[] ind1, int[] ind2)
		{
			int i;
			for(i=0; i<ind1.Length; i++)
				if(!this.PixelIncluded[ind1[i],position])
					return false;
			for(i=0; i<ind2.Length; i++)
				if(!this.PixelIncluded[ind2[i],position])
					return false;
			return true;
		}
		#region Drawing / Painting
		/// <summary>
		/// A class for drawing ribbon bounds.
		/// </summary>
		public class RibbonDrawer
		{
			#region Natural <--> View Transforms
			/// <summary>
			/// The Viewer represents the zoom window in natural coordinates.
			/// This is usually provided by the PbViewer form object in the ObjectParser namespace.
			/// </summary>
			public RectangleF ViewerRect
			{
				get { return viewerRect; }
				set
				{
					viewerRect=value;
					CalcViewerRects();
				}
			}
			protected RectangleF viewerRect=RectangleF.Empty;
			/// <summary>
			/// The natural rectangle for a ribbon object.
			/// This tightly bounds the object's graphics in natural coordinates.
			/// </summary>
			public RectangleF NaturalRect
			{
				get { return naturalRect; }
				set
				{
					naturalRect=value;
					CalcNaturalRects();
					CalcViewerRects();
				}
			}
			protected RectangleF naturalRect=RectangleF.Empty;
			public float NaturalHeight
			{
				get
				{
					return this.PlotHeight + 4.0f*this.Rib.RibbonHeight + 4.0f*this.TweenHeight;
				}
			}
			public float NaturalWidth
			{
				get
				{
					return Rib.RibbonWidth;
				}
			}
			private float left, top, width, height,l,r,t,b, whRatio;
			/// <summary>
			/// Transforms a point from the viewerRect (painting area) to natural coordinates.
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float ViewerToNatX(float x)
			{
				return (x-viewerRect.Left) * naturalRect.Width / (float)viewerRect.Width + naturalRect.Left;
			}
			/// <summary>
			/// Transforms a point from the viewerRect (painting area) to natural coordinates.
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float ViewerToNatY(float y)
			{
				return (y-viewerRect.Top) * naturalRect.Height / (float)viewerRect.Height + naturalRect.Top;
			}
			public PointF ViewerToNat(PointF pt)
			{
				return new PointF(ViewerToNatX(pt.X), ViewerToNatY(pt.Y));
			}
			public RectangleF ViewerToNat(Rectangle rect)
			{
				l = ViewerToNatX((float)rect.Left);
				r = ViewerToNatX((float)rect.Right);
				t = ViewerToNatY((float)rect.Top);
				b = ViewerToNatY((float)rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			public RectangleF ViewerToNat(RectangleF rect)
			{
				l = ViewerToNatX(rect.Left);
				r = ViewerToNatX(rect.Right);
				t = ViewerToNatY(rect.Top);
				b = ViewerToNatY(rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			/// <summary>
			/// Transforms a point from natural coordinates to the viewerRect (painting area).
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float NatToViewerX(float x)
			{
				return (x-naturalRect.Left) * (float)viewerRect.Width / naturalRect.Width + viewerRect.Left;
			}
			/// <summary>
			/// Transforms a point from natural coordinates to the viewerRect (painting area).
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float NatToViewerY(float y)
			{
				return (y-naturalRect.Top) * (float)viewerRect.Height / naturalRect.Height + viewerRect.Top;
			}
			public PointF NatToViewer(PointF pt)
			{
				return new PointF(NatToViewerX(pt.X), NatToViewerY(pt.Y));
			}
			public RectangleF NatToViewerF(RectangleF rect)
			{

				l = NatToViewerX(rect.Left);
				r = NatToViewerX(rect.Right);
				t = NatToViewerY(rect.Top);
				b = NatToViewerY(rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			public Rectangle NatToViewer(RectangleF rect)
			{
				l = NatToViewerX(rect.Left);
				r = NatToViewerX(rect.Right);
				t = NatToViewerY(rect.Top);
				b = NatToViewerY(rect.Bottom);
				return new Rectangle((int)(l+0.5f),(int)(t+0.5f),(int)(r-l+0.5f),(int)(b-t+0.5f));
			}
			public void CalcNaturalRects()
			{
				//	WARNING:  IF YOU CHANGE THIS CODE, REMEMBER TO CHANGE the NaturalHeight PROPERTY.

				if(naturalRect!=RectangleF.Empty)
				{
					top=0.0f;
					left = naturalRect.Left;

					//	PlotRect
					width = Rib.RibbonWidth;
					height = PlotHeight;
					whRatio = width/height;
					PlotRect = new RectangleF(left,top,width,height);
				
					//	TweenRects
					if(TweenRects==null)
						TweenRects = new RectangleF[4];
					//	TweenRects[0]
					top+=height;	//	below the PlotRect
					height=TweenHeight;
					TweenRects[0] = new RectangleF(left,top,width,height);

					//	RibbonRect
					top+=height;	//	below the TweenRect[0]
					height=Rib.RibbonHeight;
					RibbonRect = new RectangleF(left,top,width,height);

					//	TweenRects[1]
					top+=height;	//	below the RibbonRect
					height=TweenHeight;
					TweenRects[1] = new RectangleF(left,top,width,height);

					//	LumRect
					top+=height;	//	below the TweenRect[1]
					height=Rib.RibbonHeight;
					LumRect = new RectangleF(left,top,width,height);

					//	TweenRects[2]
					top+=height;	//	below the LumRect
					height=TweenHeight;
					TweenRects[2] = new RectangleF(left,top,width,height);

					//	YbRect
					top+=height;	//	below the TweenRect[2]
					height=Rib.RibbonHeight;
					YbRect = new RectangleF(left,top,width,height);

					//	TweenRects[3]
					top+=height;	//	below the YbRect
					height=TweenHeight;
					TweenRects[3] = new RectangleF(left,top,width,height);

					//	LumRect
					top+=height;	//	below the TweenRect[3]
					height=Rib.RibbonHeight;
					RgRect = new RectangleF(left,top,width,height);
				}

				//	WARNING:  IF YOU CHANGE THIS CODE, REMEMBER TO CHANGE the NaturalHeight PROPERTY.
			}
			public void CalcViewerRects()
			{
				if(viewerRect!=RectangleF.Empty && viewerRect.Width>0.0f && viewerRect.Height>0.0f)
				{
					PlotRectViewer = NatToViewer(PlotRect);

					int i;
					TweenRectsViewer = new RectangleF[TweenRects.Length];
					for(i=0; i<TweenRects.Length; i++)
						TweenRectsViewer[i] = NatToViewer(TweenRects[i]);

					RibbonRectViewer = NatToViewer(RibbonRect);
					LumRectViewer = NatToViewer(LumRect);
					YbRectViewer = NatToViewer(YbRect);
					RgRectViewer = NatToViewer(RgRect);
				}
			}
			#endregion
			/// <summary>
			/// The height of the blank spot between sections.
			/// </summary>
			public float TweenHeight;
			/// <summary>
			/// The height of the plot window.
			/// </summary>
			public float PlotHeight;
			/// <summary>
			/// The ribbon which this object serves.
			/// </summary>
			public Ribbon Rib;
			/// <summary>
			/// Constructs.  Not useful until NaturalRect and ViewerRect are set.
			/// </summary>
			/// <param name="Rib"></param>
			/// <param name="PlotHeight"></param>
			/// <param name="TweenHeight"></param>
			public RibbonDrawer( Ribbon Rib, float PlotHeight, float TweenHeight )
			{
				this.Rib = Rib;
				this.PlotHeight = PlotHeight;
				this.TweenHeight = TweenHeight;

				this.GridPen = new Pen(Color.Black, 1.0f);
				this.GridPen.DashPattern = new float[] {4.0f, 4.0f};
				this.GridPen.DashStyle = DashStyle.Custom;

				this.FirstPen = new Pen(Color.DodgerBlue,1.0f);
				this.SecondPen = new Pen(Color.FromArgb(255,0,0),1.0f);
				this.ScorePen = new Pen(Color.Black,1.0f);
			}
			#region Plotting
			/// <summary>
			/// The rectangle inside which plotting occurs.
			/// </summary>
			public RectangleF PlotRect, PlotRectViewer;
			public Pen FirstPen, SecondPen, GridPen, ScorePen;
			private PointF[] _pts;
			/// <summary>
			/// 
			/// </summary>
			/// <param name="g"></param>
			/// <param name="drawPspace">drawPspace  the points are specified in a p-space, where all values are
			/// within 0 to 1.  the y-axis has already been inverted to graphics space.</param>
			/// <param name="drawPen"></param>
			public void DrawPlot(object sender, System.Windows.Forms.PaintEventArgs e, PointF[] drawPspace, Pen drawPen)
			{
				if(drawPspace!=null && drawPspace.Length>1)
				{
					int i;
					float w,h,l,t;
					w=PlotRectViewer.Width-2.0f;
					h=PlotRectViewer.Height-2.0f;
					l=PlotRectViewer.Left+1.0f;
					t=PlotRectViewer.Top+1.0f;
					if(_pts==null || _pts.Length != drawPspace.Length)
						_pts= new PointF[drawPspace.Length];
					for(i=0; i<_pts.Length; i++)
						_pts[i] = new PointF( drawPspace[i].X*w+l, drawPspace[i].Y*h+t);
					e.Graphics.DrawLines(drawPen, _pts);
				}
			}
			private static float _temp;
			public Font PlotFont = new Font(FontFamily.GenericSansSerif, 10.0f);
			/// <summary>
			/// 
			/// </summary>
			/// <param name="g"></param>
			/// <param name="drawRect">This should be equal to PlotRect transformed into graphics space.</param>
			public void DrawPlot(object sender, System.Windows.Forms.PaintEventArgs e)
			{
				//	Fill with white
				e.Graphics.FillRectangle(Brushes.White, PlotRectViewer);

				//	Draw the top dash
				_temp = 0.25f*PlotRectViewer.Height;
				e.Graphics.DrawLine(this.GridPen, PlotRectViewer.Left, _temp+PlotRectViewer.Top,
					PlotRectViewer.Right, _temp+PlotRectViewer.Top);
				e.Graphics.DrawString("+pi", this.PlotFont, Brushes.Black, PlotRectViewer.Left, _temp+PlotRectViewer.Top-16.0f);

				//	Draw the middle dash
				_temp = 0.50f*PlotRectViewer.Height;
				e.Graphics.DrawLine(this.GridPen, PlotRectViewer.Left, _temp+PlotRectViewer.Top,
					PlotRectViewer.Right, _temp+PlotRectViewer.Top);

				//	Draw the lower dash
				_temp = 0.75f*PlotRectViewer.Height;
				e.Graphics.DrawLine(this.GridPen, PlotRectViewer.Left, _temp+PlotRectViewer.Top,
					PlotRectViewer.Right, _temp+PlotRectViewer.Top);
				e.Graphics.DrawString("-pi", this.PlotFont, Brushes.Black, PlotRectViewer.Left, _temp+PlotRectViewer.Top+2.0f);

				//	Draw the plots
				if(Rib.FirstPlots!=null)
				{
					int i;
					for(i=0; i<Rib.FirstPlots.Length; i++)
						DrawPlot(sender, e, Rib.FirstPlots[i], this.FirstPen);
					DrawPlot(sender, e, Rib.SecondPlot, this.SecondPen);
				}

				if(Rib.SecondPlot!=null)
				{
					//	Score the rectangle.
					e.Graphics.DrawRectangle(ScorePen, PlotRectViewer.Left, PlotRectViewer.Top,
						PlotRectViewer.Width, PlotRectViewer.Height);
				}
			}
			#endregion
			#region Drawing Tween
			public RectangleF[] TweenRects,TweenRectsViewer;
			public Brush TweenBrush = new SolidBrush(Color.Black);
			public void FillTweenRects(object sender, System.Windows.Forms.PaintEventArgs e)
			{
				for(int i=0; i<TweenRectsViewer.Length; i++)
					e.Graphics.FillRectangle(TweenBrush, TweenRectsViewer[i]);
			}
			#endregion
			#region Drawing the main ribbon.
			public RectangleF RibbonRect;
			public Rectangle RibbonRectViewer;
			public Bitmap RibbonBmp;
			public void DrawRibbon(object sender, System.Windows.Forms.PaintEventArgs e)
			{
				//	Draw the primary ribbon.
				if(RibbonBmp!=null)
					e.Graphics.DrawImage( RibbonBmp, this.RibbonRectViewer,
						0, 0, RibbonBmp.Width, RibbonBmp.Height, GraphicsUnit.Pixel );
			}
			#endregion
			#region Drawing the color ribbons.
			public RectangleF LumRect, YbRect, RgRect;
			public Rectangle LumRectViewer, YbRectViewer, RgRectViewer;
			public Bitmap LumBmp, YbBmp, RgBmp;
			public void DrawColorRibbons(object sender, PaintEventArgs e)
			{
				//	Draw the luminance ribbon.
				if(LumBmp!=null)
					e.Graphics.DrawImage( LumBmp, this.LumRectViewer, 
						0, 0, LumBmp.Width, LumBmp.Height, GraphicsUnit.Pixel );
				//	Draw the Yb ribbon.
				if(YbBmp!=null)
					e.Graphics.DrawImage( YbBmp, this.YbRectViewer, 
						0, 0, YbBmp.Width, YbBmp.Height, GraphicsUnit.Pixel );
				//	Draw the Rg ribbon.
				if(RgBmp!=null)
					e.Graphics.DrawImage( RgBmp, this.RgRectViewer, 
						0, 0, RgBmp.Width, RgBmp.Height, GraphicsUnit.Pixel );
			}
			#endregion
			public void Paint(object sender, PaintEventArgs e)
			{
				if(this.Rib.RibbonWidth>0)
				{
					this.FillTweenRects(sender, e);
					this.DrawPlot(sender, e);
					this.DrawRibbon(sender, e);
					this.DrawColorRibbons(sender, e);
				}
			}
			public static Font TextFont = new Font("Arial", 7.0f);
			public static StringFormat StrFormat;
			static RibbonDrawer()
			{
				StrFormat = new StringFormat();
				StrFormat.Alignment=StringAlignment.Far;
				StrFormat.LineAlignment=StringAlignment.Near;
				StrFormat.Trimming = StringTrimming.None;
			}
			public void DrawText(object sender, PaintEventArgs e, PointF topRight, float heightFactor)
			{
				//	Orientation and curvature
				e.Graphics.DrawString("Orientation (blue)\nCurvature (red)", TextFont, Brushes.Black, 
					topRight, StrFormat);
				//	Luminance
				topRight.Y += (LumRect.Top-PlotRect.Top)*heightFactor;
				e.Graphics.DrawString("Lum", TextFont, Brushes.Black, 
					topRight, StrFormat);
				//	Yb
				topRight.Y += (YbRect.Top-LumRect.Top)*heightFactor;
				e.Graphics.DrawString("Blue (neg)\nYellow (pos)", TextFont, Brushes.Black,
					topRight, StrFormat);
				//	Rg
				topRight.Y += (RgRect.Top-YbRect.Top)*heightFactor;
				e.Graphics.DrawString("Green (neg)\nRed (pos)", TextFont, Brushes.Black,
					topRight, StrFormat);
			}
		}
		/// <summary>
		/// The height of the main ribbon (not the total height of all the
		/// strands together, this is just the height of the main strand).
		/// </summary>
		public float RibbonHeight { get { return DirectionalDistances.Length; } }
		/// <summary>
		/// The height of a full ribbon.
		/// </summary>
		public float RibbonWidth { get { return Spine.Length; } }
		/// <summary>
		/// Paints the ribbon.
		/// </summary>
		public RibbonDrawer Painter;
		/// <summary>
		/// Calls Painter.Paint(sender, e).
		/// </summary>
		/// <param name="sender"></param>
		/// <param name="e"></param>
		public void Paint(object sender, PaintEventArgs e)
		{
			if(Painter!=null)
				Painter.Paint(sender, e);
		}
		/// <summary>
		/// Initializes the object to be painted onto a Graphics object.
		/// </summary>
		/// <param name="displayImage">The full display image.</param>
		public void InitPainting( Bitmap displayImage )
		{
			//int hlim = displayImage.Height-1;
			//int wlim = displayImage.Width-1;
			int i,j,ix,iy,ii,jj;
			float dx, dy, dxx, dyy;
			PointF pt;
			//	The height of each ribbon strand
			int ribbonHeight = DirectionalDistances.Length;
			//	Initialize the drawing object with bitmaps.
			Painter = new RibbonDrawer(this, 100.0f, 2.0f);
			if(this.RibbonWidth>1)
			{
				Painter.RibbonBmp = new Bitmap((int)(this.RibbonWidth+0.5f), (int)(this.RibbonHeight+0.5f));
				//	For creating RibbonBmp.
				float[] argb = new float[4];
				Color col;
				//	for each point along the Spine
				for(i=0; i<Spine.Length; i++)
				{
					//	for each sample along the perpendicular
					for(j=0; j<ribbonHeight; j++)
					{
						if(PixelIncluded[ribbonHeight-1-j,i])
						{
							//	The point on the image.
							pt = Pixels[ribbonHeight-1-j,i];

							//	2d interpolation
							ix = (int)pt.X;
							iy = (int)pt.Y;
							dx = pt.X-(float)ix;
							dy = pt.Y-(float)iy;
							dxx=1.0f-dx;
							dyy=1.0f-dy;

							argb[0]=argb[1]=argb[2]=argb[3]=0.0f;
							ii = ix;
							jj = iy;
							col = displayImage.GetPixel(ii,jj);
							argb[0]+=(float)col.A*dxx*dyy;
							argb[1]+=(float)col.R*dxx*dyy;
							argb[2]+=(float)col.G*dxx*dyy;
							argb[3]+=(float)col.B*dxx*dyy;

							jj = iy+1;
							col = displayImage.GetPixel(ii,jj);
							argb[0]+=(float)col.A*dxx*dy;
							argb[1]+=(float)col.R*dxx*dy;
							argb[2]+=(float)col.G*dxx*dy;
							argb[3]+=(float)col.B*dxx*dy;

							ii = ix+1;
							jj = iy;
							col = displayImage.GetPixel(ii,jj);
							argb[0]+=(float)col.A*dx*dyy;
							argb[1]+=(float)col.R*dx*dyy;
							argb[2]+=(float)col.G*dx*dyy;
							argb[3]+=(float)col.B*dx*dyy;

							jj = iy+1;
							col = displayImage.GetPixel(ii,jj);
							argb[0]+=(float)col.A*dx*dy;
							argb[1]+=(float)col.R*dx*dy;
							argb[2]+=(float)col.G*dx*dy;
							argb[3]+=(float)col.B*dx*dy;

							col = Color.FromArgb(
								(int)(argb[0]+0.5f),
								(int)(argb[1]+0.5f),
								(int)(argb[2]+0.5f),
								(int)(argb[3]+0.5f) );
							Painter.RibbonBmp.SetPixel(i,ribbonHeight-1-j,col);
						}
						else
							Painter.RibbonBmp.SetPixel(i,ribbonHeight-1-j,Color.White);
					}
				}
				//	Prepare to plot curves
				calcZeroCrossingsAndCurvePlots();

				//	Get ribbon dimensions.
				int rows=PixelValues.GetLength(0);
				int cols=PixelValues.GetLength(1);

				//	Prepare other bitmaps
				Painter.LumBmp = new Bitmap(cols,rows);
				Painter.YbBmp = new Bitmap(cols,rows);
				Painter.RgBmp = new Bitmap(cols,rows);

				//	Fill in the bitmaps using the colormap:  Util.BlueGreenGrayYellowRed
				for(i=0; i<rows; i++)
				{
					for(j=0; j<cols; j++)
					{
						if(this.PixelIncluded[i,j])
						{
							Painter.LumBmp.SetPixel(j,i,Util.BlueGreenGrayYellowRed( (LogLum[i,j]-MeanLogLum)/SdLogLum ));
							Painter.YbBmp.SetPixel(j,i,Util.BlueGreenGrayYellowRed( (LogYb[i,j]-MeanLogYb)/SdLogYb ));
							Painter.RgBmp.SetPixel(j,i,Util.BlueGreenGrayYellowRed( (LogRg[i,j]-MeanLogRg)/SdLogRg ));
						}
						else
						{
							Painter.LumBmp.SetPixel(j,i,Color.Black);
							Painter.YbBmp.SetPixel(j,i,Color.Black);
							Painter.RgBmp.SetPixel(j,i,Color.Black);
						}
					}
				}
			}
		}
		/// <summary>
		/// Indices of zero-crossing of the first angular derivative.
		/// Non-positive to positive crossings and visa versa.
		/// </summary>
		public int[] Cross1;
		/// <summary>
		/// Gathers information in CalcZeroCrossingsandPlots
		/// </summary>
		private static ArrayList z1 = new ArrayList();
		/// <summary>
		/// Plots of orientation in all visible bins of 2pi.
		/// </summary>
		public PointF[][] FirstPlots;
		/// <summary>
		/// Intraplot Coordinates of the second angular derivate plot [within (0,1; 0,1) already y-inverted]
		/// 0 is the lower plot (for wrapping values above 0), 2 is the upper plot for all other values.
		/// </summary>
		public PointF[] SecondPlot;
		/// <summary>
		/// Computes Zero1 and Zero2 from FirstAngularDerivative and SecondAngularDerivative.
		/// </summary>
		public void calcZeroCrossingsAndCurvePlots()
		{
			if(FirstAngularDerivative.Length>1)
			{
				//	Find the minimum and maximum wrapping values;
				float upCut = (float)(0.5*Math.PI);
				float downCut = -upCut;
				int wmin=0, wmax=0, wrap=0; 
				float last=(float)FirstAngularDerivative[0];
				int i;
				for(i=1; i<FirstAngularDerivative.Length; i++)
				{
					//	Downward wrap from low negative angle to high positive angle
					if(last<downCut && FirstAngularDerivative[i]>upCut)
					{
						if(--wrap < wmin)
							wmin=wrap;
					}
					//	Upward wrap from high positive angle to low negative angle
					else if(last>upCut && FirstAngularDerivative[i]<downCut)
					{
						if(++wrap > wmax)
							wmax=wrap;
					}
					last = (float)FirstAngularDerivative[i];
				}

				//	Start with one plot in each direction and one center plot.
				//	Now wmin represents the number of additional above-center plots
				//	Now wmax represents the number of additional below-center plots
				wrap = wmax-wmin+3;
				FirstPlots = new PointF[wrap][];
				float offset, fac = 1.0f/(float)(FirstAngularDerivative.Length-1);
				int j;
				for(i=0; i<wrap; i++)
				{
					FirstPlots[i]=new PointF[FirstAngularDerivative.Length];
					offset = 2.0f*(float)(i+wmin-1)*(float)Math.PI;
					for(j=0; j<FirstAngularDerivative.Length; j++)
					{
						if( j>0 )
						{
							if( FirstAngularDerivative[j-1]<downCut && FirstAngularDerivative[j]>upCut )
								offset -= (float)Util.TwoPi;
							else if( FirstAngularDerivative[j-1]>upCut && FirstAngularDerivative[j]<downCut )
								offset += (float)Util.TwoPi;
						}
						FirstPlots[i][j] = new PointF( (float)j*fac,
							Math.Min(1.0f,Math.Max(0.0f,(float)(0.5-(offset+FirstAngularDerivative[j])/4.0/Math.PI))) );
					}
				}

				//	Plot of second derivative
				SecondPlot = new PointF[SecondAngularDerivative.Length];
				for(i=0; i<SecondAngularDerivative.Length; i++)
					SecondPlot[i] = new PointF( (float)i*fac,
						Math.Min(1.0f,Math.Max(0.0f,0.5f-(float)(SecondAngularDerivative[i]/Math.PI))) );


				/*
				z1.Clear();
				bool pos1, pos1n;
				int i=0, ii, jj, k, kk;
				pos1 = FirstAngularDerivative[i]>0.0f;
				//	Count the number of respective zero crossings.
				for(i=1; i<FirstAngularDerivative.Length; i++)
				{
					pos1n = FirstAngularDerivative[i]>0.0f;
					if(pos1n!=pos1)
					{
						z1.Add(i);
						pos1=pos1n;
					}
				}

				//	Determines a cross.
				Cross1 = (int[])z1.ToArray(typeof(int));

				//	Initialize plot memory
				//	Plot of first derivative
				FirstPlot = new PointF[3][][];
				FirstPlot[0] = new PointF[Cross1.Length+1][];
				FirstPlot[1] = new PointF[Cross1.Length+1][];
				FirstPlot[2] = new PointF[Cross1.Length+1][];

				//	Plot of second derivative
				SecondPlot = new PointF[SecondAngularDerivative.Length];

				float fac = 1.0f/(float)(FirstAngularDerivative.Length-1);
				PointF pt;
				bool down;

				jj=0;
				//	For each crossing (between "bins"), plus the last "bin"
				//	It's not really a bin, but it's a plotting point, and crossings separate these points.
				for(i=0; i<=Cross1.Length; i++)
				{
					ii=jj;
					if(i==Cross1.Length)
						jj = FirstAngularDerivative.Length-1;
					else
						jj=Cross1[i]-1;

					down = FirstAngularDerivative[ii]>0.0;

					//	Fill in the middle plot
					FirstPlot[1][i] = new PointF[jj-ii+1];
					//	Fill either the upper or lower plot
					if(down) //	lower
						FirstPlot[0][i] = new PointF[jj-ii+1];
					else //	upper
						FirstPlot[2][i] = new PointF[jj-ii+1];

					for(k=ii; k<=jj; k++)
					{
						kk = k-ii;
						//	Render x and y in 0 to 1, invert for screen
						//pt = FirstPlot[1][i][kk] = new PointF((float)k*fac,
						//	(float)(0.5-FirstAngularDerivative[i]/4.0/Math.PI) );
						pt = FirstPlot[1][i][kk] = new PointF((float)k*fac,
							(float)(0.5-FirstAngularDerivative[kk]/Math.PI/2.0));
						if(down)
						{
							//	shift down (up in screen coordinates, intraplot)
							pt.Y+=0.5f;
							FirstPlot[0][i][kk] = pt;
						}
						else
						{
							//	shift up (down in screen coordinates, intraplot)
							pt.Y-=0.25f;
							FirstPlot[2][i][kk] = pt;
						}
					}
				}
				jj=0;
				//	For each crossing (between "bins"), plus the last "bin"
				//	It's not really a bin, but it's a plotting point, and crossings separate these points.
				for(i=0; i<SecondAngularDerivative.Length; i++)
					SecondPlot[i] = new PointF((float)i*fac,
						Math.Min(1.0f,Math.Max(0.0f,0.5f-(float)SecondAngularDerivative[i])));
				*/
			}
		}
		#endregion
	}
	/// <summary>
	/// Represents a single visible surface for a single parsed object.
	/// </summary>
	public class SurfaceBound
	{
		/// <summary>
		/// The ParsedObject on which the visible surface sits.
		/// </summary>
		public ParsedObject Po;
		/// <summary>
		/// The ordered elements that define the boundary of the visible surface.  This contains
		/// only OcclusionBridge objects or Edge objects in an alternating order (usually starts with an OcclusionBridge).
		/// </summary>
		public object[] Elements;
		/// <summary>
		/// The spine of each ribbon which bounds this surface.  This is
		/// computed in the constructor.
		/// </summary>
		public PointF[][] Spines;
		/// <summary>
		/// The PathF object associated with each spine.
		/// </summary>
		public PathF[] PfSpines;
		/// <summary>
		/// The ribbons which bound this surface.
		/// This is computed when CalcRibbons(...) is called.
		/// </summary>
		public Ribbon[] Ribbons;
		/// <summary>
		/// The occlusion bridges associated with this surface.  This is computed
		/// in the constructor.
		/// </summary>
		public OcclusionBridge[] Bridges;
		/// <summary>
		/// The number of occlusion bridges associated with this surface.
		/// </summary>
		public int Nbridges
		{
			get
			{
				if(Bridges==null)
					return 0;
				return Bridges.Length;
			}
		}
		/// <summary>
		/// The precision of the spines in pixels.  Consecutive elements of the spine are
		/// separated by this distance apart.
		/// </summary>
		public float Dpix;
		/// <summary>
		/// The core pixels in the entire image (integers).  This is set from the call to CalcRibbions(...).
		/// </summary>
		public Ushort3 CorePix;
		/// <summary>
		/// A factor for converting elements of CorePix to an absolute RGB value.
		/// </summary>
		public float RgbFac;
		/// <summary>
		/// The directional distances associated with each ribbon.  This is set from the call
		/// to CalcRibbons(...).
		/// </summary>
		public float[] DirDistances;
		/// <summary>
		/// The standard deviation of the curvature smoothing kernel used for constructing
		/// ribbons. This is set by the call to CalcRibbons(...).
		/// </summary>
		public float SdCurvatureKernel;
		/// <summary>
		/// Determines whether each element goes forward.  Created in constructor (alongside Elements).
		/// </summary>
		public bool[] ElementGoesForward;
		/// <summary>
		/// Tells if the surface is valid (or not).  This value is set in the constructor.
		/// </summary>
		public bool IsValidSurface
		{
			get{ return isValidSurface; }
		}
		/// <summary>
		/// Determined in constructor.
		/// </summary>
		private bool isValidSurface=false;
		/// <summary>
		/// Constructs a SurfaceBound of the ParsedObject given a set of 
		/// ownedEdges, bridges, and unownedEdges which are not already part of another SurfaceBound.
		/// Most objects will only have one SurfaceBound unless they are self-occluding.
		/// On exit, ownedEdges, bridges, and unownedEdges will contain information about any remaining
		/// surfaces (to construct an additional SurfaceBound if necessary).
		/// Specifically, if (on exit) ownedEdges has a count > 0, then the object is self-occluding.
		/// </summary>
		/// <param name="Po"></param>
		/// <param name="ownedEdges"></param>
		/// <param name="bridges"></param>
		/// <param name="unownedEdges"></param>
		public SurfaceBound( ParsedObject Po, float Dpix,
			EdgeArray ownedEdges, OcclusionBridgeList bridges, EdgeArray unownedEdges )
		{
			this.Po = Po;
			this.Dpix = Dpix;
			ArrayList els = new ArrayList();

			if(ownedEdges==null || ownedEdges.Count==0)
				isValidSurface=false;
			else
			{
				Edge ed;
				OcclusionBridge br, br2;
				VertexF vfSub,vfClose;
				int i;
				//	The direction of the path with respect to the edges on the polygon.
				//	If true, move from Begin-->End.  Otherwise, move End-->Begin.
				//	This value must be determined after a bridge is crossed, but for starters,
				//	we will grow forward.
				bool growForward = true;
				//	The ownedEdge that starts us off at the point of subjection.
				//	Herein, the 'subject' is always the point where the surface boundary is growing.
				object subject = ownedEdges[ownedEdges.Count-1];
				ownedEdges.RemoveAt(ownedEdges.Count-1);
				els.Add(subject);
				vfClose = ((Edge)subject).Begin;	//	The closure vertex
				vfSub = ((Edge)subject).End;			//	The vertex at the point of subjection
				//	The polygon at the point of subjection
				//	(i.e. the next element must be attached to this polygon).
				PolygonF pfSub = ((Edge)subject).Owner;
				//	Says that an edge has been found.
				bool found;
				bool fromUnowned=false;
				//	Determines when we're done growing
				bool done = false;
				while(!done)
				{
					found=false;
					//	For each owned edge
					for(i=0; i<ownedEdges.Count && !found; i++)
					{
						ed = ownedEdges[i];
						//	Where the edge is part of pfSub
						if(pfSub==ed.Owner)
						{
							//	And where we are growing from an owned edge
							if(subject.GetType()==typeof(Edge))
							{
								//	And if a joining candidate exists
								if( (growForward&&(vfSub==ed.Begin)) || ((!growForward)&&(vfSub==ed.End)) )
								{
									//	A new subject has been found.
									found=true;
									subject=ed;
									els.Add(ed);
									ownedEdges.RemoveAt(i);
									//	Set the new point of subjection.
									if(growForward)
										vfSub=ed.End;
									else
										vfSub=ed.Begin;
									//	For next round, we are not growing from an unowned edge.
									fromUnowned=false;
								}
							}
								//	And where we are growing from an occlusion bridge
							else
							{
								//	And if a joining candidate exists
								if(vfSub==ed.Begin||vfSub==ed.End)
								{
									//	TO DO:  The best joining candidate is an edge from the same polygon.
									//	If joining to a different polygon, make sure that the better candidate
									//	does not exist in the unsearched portion of ownedEdges.

									//	A new subject has been found.
									found=true;
									subject=ed;
									els.Add(ed);
									ownedEdges.RemoveAt(i);
									//	Determine directio of growth
									growForward = vfSub==ed.Begin;
									//	Set the new point of subjection.
									if(growForward)
										vfSub=ed.End;
									else
										vfSub=ed.Begin;
									//	For next round, we are not growing from an unowned edge.
									fromUnowned=false;
								}
							}
						}
					}
					//	If a joining edge was not found, and if we are growing from
					//	an owned or unowned edge, then look for an occlusion bridge to grow onto.
					if(!found && (subject.GetType()==typeof(Edge) || fromUnowned))
					{
						//	For each bridge
						for(i=0; i<bridges.Count && !found; i++)
						{
							br = bridges[i];
							//	The second element is a joining candidate, so reverse the bridge.
							//	By making sure the bridge is reversed, we guarantee that els
							//	contains all bridges going in the correct order around the polygon.
							if(br.Pf2==pfSub && br.Vf2==vfSub)
							{
								br.Pf2=br.Pf1;
								br.Vf2=br.Vf1;
								br.Pf1=pfSub;
								br.Vf1=vfSub;
								//	The next "if" statement will be followed through.
							}
							//	If the first element is a joining candidate
							if(br.Pf1==pfSub && br.Vf1==vfSub)
							{
								//	A new subject has been found.
								found=true;
								subject=br;
								els.Add(br);
								bridges.RemoveAt(i);
								//	Set the new point of subjection.
								vfSub=br.Vf2;
								pfSub=br.Pf2;
								//	For next round, we are not growing from an unowned edge.
								fromUnowned=false;
							}
						}
					}
					//	If a joining section was still not found, and if we are growing from an owned edge,
					//	then try to grow along unowned edges (and define a new bridge by that path).
					if(!found && subject.GetType()==typeof(Edge))
					{
						//	The latent occlusion bridge.
						br = new OcclusionBridge(pfSub, vfSub, pfSub, null);
						while(!found)
						{
							//	We are finished growing unless an unowned edge is found.
							found=true;
							//	Try to find find an unowned edge bound by vfSub
							for(i=0; i<unownedEdges.Count; i++)
							{
								ed=unownedEdges[i];
								//	And if a joining candidate exists
								if( (growForward&&(vfSub==ed.Begin)) || ((!growForward)&&(vfSub==ed.End)) )
								{
									//	Remove the unowned edge from future consideration
									unownedEdges.RemoveAt(i);
									//	Set the new point of subjection.
									if(growForward)
										vfSub=ed.End;
									else
										vfSub=ed.Begin;
									//	Not finished growing
									found=false;
									//	Next time will be from an unowned edge
									fromUnowned=true;
								}
							}
						}
						//	If vfSub has changed, then update the new occlusion bridge and add to solution
						if(vfSub!=br.Vf1)
						{
							br.Vf2=vfSub;

							//	A new subject has been found.
							found=true;
							subject=br;
							els.Add(br);
							//	Set the new point of subjection.
							vfSub=br.Vf2;
							//	For next round, we are not growing from an unowned edge.
							fromUnowned=true;
						}
							//	Otherwise, growing has completely stopped.
						else
							done=true;
					}
					//	Finished growing
					if(!found || vfSub==vfClose)
						done=true;
				}
				object ob1, ob2=els[0];
				//	Growing has finished.  Join any consecutive occlusion bridges.
				for(i=1; i<els.Count; i++)
				{
					ob1=ob2;
					ob2 = els[i];
					if(ob1.GetType()==typeof(OcclusionBridge) && ob2.GetType()==typeof(OcclusionBridge))
					{
						br = (OcclusionBridge)ob1;
						br2= (OcclusionBridge)ob2;
						els.RemoveAt(i);
						ob2=els[i-1]=new OcclusionBridge(br.Pf1,br.Vf1,br2.Pf2,br2.Vf2);
						i--;
					}
				}
				
				//	The possible solution is contained in els.
				//	This is only a real solution if there is closure.  There is closure
				//	if the end of the last element maps to the beginning of the first element.
				object lastEl=els[els.Count-1];
				ed = (Edge)els[0];
				//	Easy, the last element is an occlusion bridge.
				if(lastEl.GetType()==typeof(OcclusionBridge))
					isValidSurface = ((OcclusionBridge)lastEl).Vf2==ed.Begin;
				else
					isValidSurface = ((Edge)lastEl).End==ed.Begin;

				//	Order with an occlusion bridge to start (if exists), but only if this is a valid surface.
				if(isValidSurface)
				{
					//	Find the last occlusion bridge (if exists)
					int ind = -1;
					for(i=els.Count-1; i>0 && ind<0; i--)
						if(els[i].GetType()==typeof(OcclusionBridge))
							ind=i;
					//	Order the elements
					Elements = new object[els.Count];
					if(ind>0)
					{
						i=0;
						while(i<Elements.Length)
						{
							Elements[i]=els[ind++];
							if(ind>=Elements.Length)
								ind=0;
							i++;
						}
					}
					else
						Elements = els.ToArray();
				}
				else
					Elements = els.ToArray();

				//	Determine the direction of growth.  Guaranteed the first element grows forward
				//	because it is either a bridge (all of which go forward) or the first chosen edge element.
				ElementGoesForward = new bool[Elements.Length];
				ElementGoesForward[0]=true;
				for(i=1; i<Elements.Length; i++)
				{
					if(ElementGoesForward[i-1])
						ElementGoesForward[i]=EndVf(Elements[i-1])==BeginVf(Elements[i]);
					else
						ElementGoesForward[i]=EndVf(Elements[i])!=BeginVf(Elements[i-1]);
				}

				//	Calculate the spines
				calcSpinesAndBridges();
			}
		}
		/// <summary>
		/// Called in the constructor.  Constructs spines from Edge objects (found in Elements) and
		/// also compiles a list of bridges (also found in Elements).
		/// </summary>
		private void calcSpinesAndBridges()
		{
			if(Elements!=null && Elements.Length>0)
			{
				Edge ed;
				PathF pf;
				VertexF vf1, vf2;
				//	Count the number of spines
				int ct = 0, ctBridge=0;
				bool onEdge = Elements[0].GetType()==typeof(Edge);
				if(onEdge)
					ct=1;
				else
					ctBridge=1;
				for(int i=1; i<Elements.Length; i++)
				{
					if(!onEdge && Elements[i].GetType()==typeof(Edge))
					{
						onEdge = true;
						ct++;
					}
					else if(onEdge && Elements[i].GetType()==typeof(OcclusionBridge))
						onEdge = false;
					if(Elements[i].GetType()==typeof(OcclusionBridge))
						ctBridge++;
				}
				//	Allocate memory
				Spines = new PointF[ct][];
				PfSpines = new PathF[ct];
				Bridges = new OcclusionBridge[ctBridge];
				ct=-1;
				ctBridge=-1;
				float carryOver=0.0f;
				//	Get the spines
				for(int i=0; i<Elements.Length; i++)
				{
					//	If we've found the beginning of a spine
					if(Elements[i].GetType()==typeof(Edge))
					{
						ed = (Edge)Elements[i];
						pf = ed.Owner;

                        //  See if it's a full enclosure
                        if( Elements.Length==1 && ed.IsFullEnclosure )
                        {
                            vf1=vf2=pf[0];
                        }
                        else
                        {
                            //	Get the first vertex
                            if(ElementGoesForward[i])
                                vf1 = ed.Begin;
                            else
                                vf1 = ed.End;
                            while(i<Elements.Length-1 && Elements[i+1].GetType()==typeof(Edge))
                                i++;
                            ed = (Edge)Elements[i];
                            //	Get the last vertex, note, ElementGoesForward[i] will not have changed
                            if(ElementGoesForward[i])
                                vf2 = ed.End;
                            else
                                vf2 = ed.Begin;
                        }

						//	Get the spine
						Spines[++ct] = pf.EdgePoints(Dpix, vf1, vf2, ElementGoesForward[i], ref carryOver);
						PfSpines[ct] = pf;
					}
					else
					{
						Bridges[++ctBridge]=(OcclusionBridge)Elements[i];
					}
				}
			}
		}
        /// <summary>
        /// This is called in order to build ribbons.
        /// </summary>
        /// <param name="CorePixels">Linear pixel values of the image, as measured
        /// by the calibrated Kodak DCS720x Camera, in milliwatts per steridian per square meter.</param>
        /// <param name="DisplayBmp">The bitmap object used to display the image.</param>
        /// <param name="DirDistances">The height to width ratio of each side of the ribbon.</param>
        /// <param name="SdCurvatureKernel"></param>
        public void CalcRibbons(Ushort3 CorePix,float RgbFac,double heightProportion,float SdCurvatureKernel,ColorSpace[] spaces)
        {
            int j,ind;

            //	Store information
            this.CorePix = CorePix;
            this.RgbFac = RgbFac;
            
            this.SdCurvatureKernel = SdCurvatureKernel;

            int sideLen = (int)(Math.Sqrt((double)Po.Area)*heightProportion);
            DirDistances = new float[2*sideLen+1];
            DirDistances[sideLen] = 0.0f;
            for(ind=0;ind<sideLen;ind++)
            {
                DirDistances[ind] = ind-sideLen;
                DirDistances[2*sideLen-ind] = sideLen-ind;
            }

            //	Calculate ribbons
            Ribbons = new Ribbon[Spines.Length];
            for(int i=0;i<Ribbons.Length;i++)
            {
                Ribbons[i] = new Ribbon(
                    Po,PfSpines[i],CorePix,RgbFac,(double)Dpix,Spines[i],DirDistances,SdCurvatureKernel,null,0,null);
                for(j=0;j<spaces.Length;j++)
                    Ribbons[i].CalcColorRibbons(spaces[j]);
            }
        }
		/// <summary>
		/// This is called in order to build ribbons.
		/// </summary>
		/// <param name="CorePixels">Linear pixel values of the image, as measured
		/// by the calibrated Kodak DCS720x Camera, in milliwatts per steridian per square meter.</param>
		/// <param name="DisplayBmp">The bitmap object used to display the image.</param>
		/// <param name="DirDistances">The directional distance values of the ribbon.  Negative values
		/// outside polygon, positive inside, with magnitude representing distance from the edge.</param>
		/// <param name="SdCurvatureKernel"></param>
		public void CalcRibbons( Ushort3 CorePix, float RgbFac, float[] DirDistances, float SdCurvatureKernel, 
            ColorSpace[] spaces,ComplexD[][,] filters,int filterLength, PointF[] filterOffsets)
		{
			int j;

			//	Store information
			this.CorePix = CorePix;
			this.RgbFac = RgbFac;
			this.DirDistances = DirDistances;
			this.SdCurvatureKernel = SdCurvatureKernel;

			//	Calculate ribbons
			Ribbons = new Ribbon[Spines.Length];
			for(int i=0; i<Ribbons.Length; i++)
			{
				Ribbons[i] = new Ribbon(
					Po, PfSpines[i], CorePix, RgbFac, (double)Dpix, Spines[i], DirDistances, SdCurvatureKernel, filters, filterLength, filterOffsets );
				for(j=0; j<spaces.Length; j++)
					Ribbons[i].CalcColorRibbons(spaces[j]);
			}
		}
		/// <summary>
		/// Makes a single spine denoting the entire surface bound, interpolating through occlusion bridges.
		/// </summary>
		/// <param name="dpix">The distance between spine points.</param>
		/// <returns></returns>
		public PointF[] MakeSingleSpine(float dpix)
		{
			PointF[] output;
			ArrayList agg = new ArrayList();
			int len = 0;

			float carryOver = dpix;

			if(!isValidSurface)
				return null;

			int i;
			for(i=0; i<Elements.Length; i++)
			{
				if(Elements[i].GetType()==typeof(OcclusionBridge))
				{
					output = ((OcclusionBridge)Elements[i]).EdgePoints(dpix, ElementGoesForward[i], ref carryOver);
					if(output!=null && output.Length>0)
					{
						agg.Add(output);
						len+=output.Length;
					}
				}
				else if(Elements[i].GetType()==typeof(Edge))
				{
					output= ((Edge)Elements[i]).EdgePoints(dpix, ref carryOver);
					if(output!=null && output.Length>0)
					{
						agg.Add(output);
						len+=output.Length;
					}
				}
			}

			output = new PointF[len];
			PointF[] temp;
			int ct=0,j;
			for(i=0; i<agg.Count; i++)
			{
				temp = ((PointF[])agg[i]);
				for(j=0; j<temp.Length; j++)
					output[ct+j] = temp[j];
				ct+=temp.Length;
			}
			return output;
		}
		/// <summary>
		/// Used to extract a beginning VertexF from either an Edge or OcclusionBridge object.
		/// </summary>
		/// <param name="o">An Edge or OcclusionBridge object.</param>
		/// <returns></returns>
		public static VertexF BeginVf(object o)
		{
			if(o.GetType()==typeof(Edge))
				return ((Edge)o).Begin;
			else if(o.GetType()==typeof(OcclusionBridge))
				return ((OcclusionBridge)o).Vf1;
			return null;
		}
		/// <summary>
		/// Used to extract an ending VertexF from either an Edge or OcclusionBridge object.
		/// </summary>
		/// <param name="o">An Edge or OcclusionBridge object.</param>
		/// <returns></returns>
		public static VertexF EndVf(object o)
		{
			if(o.GetType()==typeof(Edge))
				return ((Edge)o).End;
			else if(o.GetType()==typeof(OcclusionBridge))
				return ((OcclusionBridge)o).Vf2;
			return null;
		}
		#region Drawing / Painting
		/// <summary>
		/// A class for drawing surface bounds.
		/// </summary>
		public class SurfaceDrawer
		{
			#region Natural <--> View Transforms
			public RectangleF ViewerRect
			{
				get
				{
					return viewerRect;
				}
				set
				{
					viewerRect=value;
					//	Calculate the view rectangles
					CalcViewRects();
					//	Pass the view rectangles to the ribbons
					PassViewRects();
				}
			}
			private RectangleF viewerRect;
			public RectangleF NaturalRect
			{
				get{ return naturalRect; }
				set
				{
					naturalRect = value;
					NaturalOffset.X=naturalRect.X;
					NaturalOffset.Y=naturalRect.Y;
					//	Calculate the natural rectangles
					CalcNaturalRects();
					//	Pass the natural rectangles to the ribbons.
					PassNaturalRects();
					//	Calculate the view rectangles
					CalcViewRects();
					//	Pass the view rectangles to the ribbons
					PassViewRects();
				}
			}
			private RectangleF naturalRect;
			public PointF NaturalOffset;
			public float SumRibbonWidths
			{
				get
				{
					if(Surf==null || Surf.Ribbons==null || Surf.Ribbons.Length==0)
						return 0.0f;
					float output=0.0f;
					for(int i=0; i<Surf.Ribbons.Length; i++)
						output+=Surf.Ribbons[i].RibbonWidth;
					return output;
				}
			}
			public float SumBridgeWidths
			{
				get
				{
					if(Surf==null || Surf.Ribbons==null || Surf.Ribbons.Length==0)
						return 0.0f;
					float output=0.0f;
					for(int i=0; i<Surf.Bridges.Length; i++)
						output+=(float)Surf.Bridges[i].BridgeLength;
					return output;
				}
			}
			public float NaturalWidth
			{
				get
				{
					if(Surf==null || Surf.Ribbons==null || Surf.Ribbons.Length==0)
						return 0.0f;
					//return (float)Surf.Nbridges*BridgeWidth+SumRibbonWidths;
					return SumBridgeWidths+SumRibbonWidths;
				}
			}
			public float NaturalHeight
			{
				get
				{
					if(Surf==null || Surf.Ribbons==null || Surf.Ribbons.Length==0)
						return 0.0f;
					return this.Surf.Ribbons[0].Painter.NaturalHeight;
				}
			}
			private float l,r,t,b;
			/// <summary>
			/// Transforms a point from the viewerRect (painting area) to natural coordinates.
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float ViewerToNatX(float x)
			{
				return (x-viewerRect.Left) * naturalRect.Width / (float)viewerRect.Width + naturalRect.Left;
			}
			/// <summary>
			/// Transforms a point from the viewerRect (painting area) to natural coordinates.
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float ViewerToNatY(float y)
			{
				return (y-viewerRect.Top) * naturalRect.Width / (float)viewerRect.Width + naturalRect.Top;
			}
			public PointF ViewerToNat(PointF pt)
			{
				return new PointF(ViewerToNatX(pt.X), ViewerToNatY(pt.Y));
			}
			public RectangleF ViewerToNat(Rectangle rect)
			{
				l = ViewerToNatX((float)rect.Left);
				r = ViewerToNatX((float)rect.Right);
				t = ViewerToNatY((float)rect.Top);
				b = ViewerToNatY((float)rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			public RectangleF ViewerToNat(RectangleF rect)
			{
				l = ViewerToNatX(rect.Left);
				r = ViewerToNatX(rect.Right);
				t = ViewerToNatY(rect.Top);
				b = ViewerToNatY(rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			/// <summary>
			/// Transforms a point from natural coordinates to the viewerRect (painting area).
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float NatToViewerX(float x)
			{
				return (x-naturalRect.Left) * (float)viewerRect.Width / naturalRect.Width + viewerRect.Left;
			}
			/// <summary>
			/// Transforms a point from natural coordinates to the viewerRect (painting area).
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float NatToViewerY(float y)
			{
				return (y-naturalRect.Top) * (float)viewerRect.Width / naturalRect.Width + viewerRect.Top;
			}
			public PointF NatToViewer(PointF pt)
			{
				return new PointF(NatToViewerX(pt.X), NatToViewerY(pt.Y));
			}	
			public RectangleF NatToViewerF(RectangleF rect)
			{
				l = NatToViewerX(rect.Left);
				r = NatToViewerX(rect.Right);
				t = NatToViewerY(rect.Top);
				b = NatToViewerY(rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			public Rectangle NatToViewer(RectangleF rect)
			{
				l = NatToViewerX(rect.Left);
				r = NatToViewerX(rect.Right);
				t = NatToViewerY(rect.Top);
				b = NatToViewerY(rect.Bottom);
				return new Rectangle((int)(l+0.5f),(int)(t+0.5f),(int)(r-l+0.5f),(int)(b-t+0.5f));
			}
			#endregion
			/// <summary>
			/// The surface bound that this object serves.
			/// </summary>
			public SurfaceBound Surf;
			/// <summary>
			/// The natural rectangles for each ribbon.
			/// </summary>
			public RectangleF[] RibbonRects, RibbonRectsViewer;
			/// <summary>
			/// The natural rectangles for each empty space between surfaces.
			/// </summary>
			public RectangleF[] BridgeRects, BridgeRectsViewer;
			/// <summary>
			/// The brush used to fill occlusion bridges.
			/// </summary>
			public Brush BridgeBrush = new SolidBrush(Color.LightGray);
			/// <summary>
			/// Initializes a SurfaceDrawer object.
			/// </summary>
			/// <param name="NaturalOffset">The offset of this surface in natural coordinates..</param>
			/// <param name="Surf">The surface bound upon which this object is to be based.</param>
			public SurfaceDrawer(PointF NaturalOffset, SurfaceBound Surf)
			{
				this.Surf=Surf;
				this.NaturalOffset=NaturalOffset;
				this.naturalRect=new RectangleF(NaturalOffset.X, NaturalOffset.Y, NaturalWidth, NaturalHeight);

				if(Surf!=null && Surf.Ribbons!=null && Surf.Ribbons.Length>0)
				{
					//	Create a surface rect for each Surfaces item.
					RibbonRects = new RectangleF[Surf.Ribbons.Length];
					BridgeRects = new RectangleF[Surf.Nbridges];

					//	Allocate memory.
					RibbonRectsViewer = new RectangleF[Surf.Ribbons.Length];
					BridgeRectsViewer = new RectangleF[Surf.Nbridges];
				}

				//	Calculate the natural rectangles
				CalcNaturalRects();
				//	Pass the natural rectangles to the ribbons.
				PassNaturalRects();
			}
			/// <summary>
			/// Calculates the natural rectangles for ribbons.
			/// </summary>
			public void CalcNaturalRects()
			{
				if(Surf!=null && Surf.Ribbons!=null && Surf.Ribbons.Length>0)
				{
					int i;
					float left=naturalRect.Left;
					//	If Elements[0] is a bridge, then do bridges first.  Otherwise, do them last.
					bool bridgeFirst = Surf.Elements[0].GetType()==typeof(OcclusionBridge);
					//	For each ribbon
					for(i=0; i<Surf.Ribbons.Length; i++)
					{
						//	Bridge first
						if(bridgeFirst && i<Surf.Bridges.Length)
						{
							BridgeRects[i] = new RectangleF(left, naturalRect.Y,
								(float)Surf.Bridges[i].BridgeLength, Surf.Ribbons[i].Painter.NaturalHeight);
							//	Increment left
							left=BridgeRects[i].Right;
						}
						//	Compute a ribbon rect
						RibbonRects[i]=new RectangleF(left, naturalRect.Y,
							Surf.Ribbons[i].RibbonWidth, Surf.Ribbons[i].Painter.NaturalHeight);
						//	Increment left
						left=RibbonRects[i].Right;
						//	Bridge after
						if(!bridgeFirst && i<Surf.Bridges.Length)
						{
							BridgeRects[i] = new RectangleF(left, naturalRect.Y,
								(float)Surf.Bridges[i].BridgeLength, Surf.Ribbons[i].Painter.NaturalHeight);
							//	Increment left
							left=BridgeRects[i].Right;
						}
					}
				}
			}
			/// <summary>
			/// Passes the natural rects for ribbons to the ribbons.
			/// CalcNaturalRects() should be called first.
			/// </summary>
			public void PassNaturalRects()
			{
				if(Surf!=null && Surf.Ribbons!=null && Surf.Ribbons.Length>0)
				{
					for(int i=0; i<RibbonRects.Length; i++)
						this.Surf.Ribbons[i].Painter.NaturalRect=RibbonRects[i];
				}
			}
			/// <summary>
			/// Calculates the natural rects for ribbons.
			/// </summary>
			public void CalcViewRects()
			{
				if(Surf!=null && Surf.Ribbons!=null && Surf.Ribbons.Length>0)
				{
					for(int i=0; i<BridgeRects.Length; i++)
						this.BridgeRectsViewer[i]=this.NatToViewerF(this.BridgeRects[i]);
					for(int i=0; i<RibbonRects.Length; i++)
						this.RibbonRectsViewer[i]=this.NatToViewerF(this.RibbonRects[i]);
				}
			}
			/// <summary>
			/// Passes the view rects for ribbons to the ribbons.
			/// CalcViewRects() should be called first.
			/// </summary>
			public void PassViewRects()
			{
				if(Surf!=null && Surf.Ribbons!=null && Surf.Ribbons.Length>0)
				{
					for(int i=0; i<RibbonRects.Length; i++)
						this.Surf.Ribbons[i].Painter.ViewerRect=RibbonRectsViewer[i];
				}
			}
			/// <summary>
			/// Paint the surface with the graphics object.  Paints bridges than ribbons.
			/// </summary>
			/// <param name="sender"></param>
			/// <param name="e"></param>
			public void Paint(object sender, PaintEventArgs e)
			{
				int i;
				//	Paint the bridges.
				if(BridgeRectsViewer!=null)
					for(i=0; i<BridgeRectsViewer.Length; i++)
						e.Graphics.FillRectangle(BridgeBrush, BridgeRectsViewer[i]);
				//	Paint the ribbons
				if(Surf!=null && Surf.Ribbons!=null)
					for(i=0; i<Surf.Ribbons.Length; i++)
						if(Surf.Ribbons[i]!=null)
							Surf.Ribbons[i].Paint(sender, e);
			}
		}
		/// <summary>
		/// The bitmap used for displaying the entire image.  This is set from the call to CalcRibbons(...).
		/// </summary>
		public Bitmap DisplayBmp;
		/// <summary>
		/// For painting surfaces.
		/// </summary>
		public SurfaceDrawer Painter;
		/// <summary>
		/// Prepares the object to be painted onto a Graphics object.  This assumes that CalcRibbons
		/// has already been called.
		/// </summary>
		/// <param name="DisplayBmp">The full display image.</param>
		public void InitPainting( Bitmap DisplayBmp )
		{
			//	Hold the display bitmap
			this.DisplayBmp = DisplayBmp;
			//	Initialize subordinate ribbons for painting
			for(int i=0; i<Ribbons.Length; i++)
				Ribbons[i].InitPainting( DisplayBmp );
			//	Initializes painting (based on subordinates) and then takes
			//	over the rectangle schema of the subordinate Ribbons.
			Painter = new SurfaceDrawer(new PointF(0.0f,0.0f),this);
		}
		#endregion
	}
	/// <summary>
	/// A list of surface bounds associated with a single ParsedObject.  Most only have a single surface, but when
	/// a surface folds over (self-occlusion), there can be more surfaces.
	/// </summary>
	public class SurfaceBoundList : CollectionBase
	{
		/// <summary>
		/// This determines whether the object was segmented correctly enough to make a valid
		/// surface representation.
		/// </summary>
		public bool IsValid
		{
			get
			{
				if(this.Count==0)
					return false;
				for(int i=0; i<this.Count; i++)
					if(!this[i].IsValidSurface)
						return false;
				return true;
			}
		}
		/// <summary>
		/// The ParsedObject to which the surfaces are associated.
		/// </summary>
		public ParsedObject Po;
		/// <summary>
		/// The precision of the Ribbon spines.
		/// </summary>
		public float Dpix;
		/// <summary>
		/// Constructs the surface bounds for an object.
		/// </summary>
		/// <param name="po">The parsed object.</param>
		/// <param name="dpix">The size of the orientation smoothing kernel.</param>
		public SurfaceBoundList(ParsedObject po, float dpix) : base()
		{
			this.Po = po;
			this.Dpix = dpix;

			PolygonF pf;
			EdgeGroup eg;
			Edge ed;
			int i,j,k;
			EdgeArray ownedEdges = new EdgeArray();
			EdgeArray unownedEdges = new EdgeArray();
			//EdgeArray surfaceSplitters = new EdgeArray();

			//	Enumerate the owned and unowned edges.
			for(i=0; i<po.Polygons.Count; i++)
			{
				pf = po.Polygons[i];
				pf.CalcEdges();
				for(j=0; j<pf.EdgeGroups.Count; j++)
				{
					eg = pf.EdgeGroups[j];
					ed = eg.DominantEdge;
					if(ed!=null && ed.Owner==pf)
					{
						ownedEdges.Add(ed);
					}
					else
					{
						for(k=0;k<eg.Count;k++)
						{
							if(k!=eg.DominantInd)
							{
								ed = eg[k];
								if(ed!=null && ed.Owner==pf)
									unownedEdges.Add(ed);
							}
						}
					}
				}
			}
			//	Enumerate the bridges.
			OcclusionBridgeList bridges = new OcclusionBridgeList();
			for(i=0; i<po.Bridges.Count; i++)
				bridges.Add(po.Bridges[i]);
			//	Construct the surface bounds
			while(ownedEdges.Count>0)
				this.Add(new SurfaceBound(Po, Dpix, ownedEdges, bridges, unownedEdges));
		}
        /// <summary>
        /// Computes the ribbons.
        /// </summary>
        /// <param name="CorePixels"></param>
        /// <param name="DirDistances"></param>
        /// <param name="SdCurvatureKernel"></param>
        public void CalcRibbons(Ushort3 corePix,float rgbFac,
            double heightProportion,float SdCurvatureKernel,ColorSpace[] spaces)
        {
            for(int i=0;i<this.Count;i++)
                this[i].CalcRibbons(corePix,rgbFac,heightProportion,SdCurvatureKernel,spaces);
        }
		/// <summary>
		/// Computes the ribbons.
		/// </summary>
		/// <param name="CorePixels"></param>
		/// <param name="DirDistances"></param>
		/// <param name="SdCurvatureKernel"></param>
		public void CalcRibbons( Ushort3 corePix, float rgbFac,
            float[] DirDistances,float SdCurvatureKernel,ColorSpace[] spaces,ComplexD[][,] filters,int filterLength,PointF[] filterOffsets)
		{
			for(int i=0; i<this.Count; i++)
				this[i].CalcRibbons(corePix, rgbFac, DirDistances, SdCurvatureKernel, spaces, filters, filterLength, filterOffsets);
		}
		/// <summary>
		/// Finds the largest section by searching the Ribbons that are part of this
		/// SurfaceBoundList.
		/// </summary>
		/// <param name="rib">The ribbon which owns the section.</param>
		/// <param name="pos">The indices demarcating the best section.</param>
		/// <returns>true if a section was found.</returns>
		public bool BestSection(out Ribbon rib, out int[] pos)
		{
			rib=null;
			pos=null;
			int maxSize=0, size;
			int i,j,k;
			int[][] currSecs;
			Ribbon currRib;

			// For each SurfaceBound
			for(i=0; i<this.Count; i++)
			{
				//	For each Ribbon on that bound
				for(j=0; j<this[i].Ribbons.Length; j++)
				{
					currRib = this[i].Ribbons[j];
					//	Get the valid sections
					currSecs=currRib.ValidSections(maxSize+1);
					//	If any are big enough
					if(currSecs.Length>0)
					{
						//	Go through all the sections
						for(k=0; k<currSecs.Length; k++)
						{
							//	Compute the size of each section
							size=currSecs[k][1]-currSecs[k][0]+1;
							//	If the size is larger than max, then this is the new max.
							if(size>maxSize)
							{
								pos=currSecs[k];
								rib=currRib;
								maxSize=size;
							}
						}
					}
				}
			}
			return maxSize>0;
		}
		/// <summary>
		/// Finds the largest section by searching the Ribbons that are part of this
		/// SurfaceBoundList.
		/// </summary>
		/// <param name="ind1"></param>
		/// <param name="ind2"></param>
		/// <param name="rib">The ribbon which owns the section.</param>
		/// <param name="pos">The indices demarcating the best section.</param>
		/// <returns>true if a section was found.</returns>
		public bool BestSection(int[] ind1, int[] ind2, out Ribbon rib, out int[] pos)
		{
			rib=null;
			pos=null;
			int maxSize=0, size;
			int i,j,k;
			int[][] currSecs;
			Ribbon currRib;

			// For each SurfaceBound
			for(i=0; i<this.Count; i++)
			{
				//	For each Ribbon on that bound
				for(j=0; j<this[i].Ribbons.Length; j++)
				{
					currRib = this[i].Ribbons[j];
					//	Get the valid sections
					currSecs=currRib.ValidSections(ind1, ind2, maxSize+1);
					//	If any are big enough
					if(currSecs.Length>0)
					{
						//	Go through all the sections
						for(k=0; k<currSecs.Length; k++)
						{
							//	Compute the size of each section
							size=currSecs[k][1]-currSecs[k][0]+1;
							//	If the size is larger than max, then this is the new max.
							if(size>maxSize)
							{
								pos=currSecs[k];
								rib=currRib;
								maxSize=size;
							}
						}
					}
				}
			}
			return maxSize>0;
		}
		#region Drawing / Painting
		public void InitPainting( Size PbSize, Bitmap DisplayBmp )
		{
			//	Initialize painting for each SurfaceBound
			for(int i=0; i<this.Count; i++)
				this[i].InitPainting(DisplayBmp);
			//	Take charge of the painting schema
			Painter = new SurfacesDrawer(PbSize, this);
		}
		/// <summary>
		/// Given a location in the PbRibbon picture-box control (in ObjectParser), associated image coordinates
		/// are returned.  This is usually the location along an object's edge, occlusion bridge, or inner contour.
		/// </summary>
		/// <param name="pbRibbonPt">PointF  Coordinates in the picture-box PbRibbon</param>
		/// <returns>Image coordinates or PointF.Empty if pbRibbonPt is not within the display
		/// area of this SurfaceBoundList.  Also returns PointF.Empty if the Painter object
		/// has not been initialized.</returns>
		public PointF GetImgCoordinates(PointF pbRibbonPt)
		{
			//	Return empty if Painter object has not been properly initialized
			if(Painter==null || Painter.Surfaces==null || Painter.Surfaces.Count<1)
				return PointF.Empty;

			//	Convert to natural coordinates
			pbRibbonPt = this.Painter.ViewerToNat(pbRibbonPt);

			//	Return empty if this SurfaceBoundList does not contain the aforementioned point.
			if( !Painter.NaturalRect.Contains(pbRibbonPt) )
				return PointF.Empty;

			//	Find which Element contains the focus
			int i, j;
			SurfaceBound sb;
			for(i=0; i<this.Count; i++)
			{
				sb = this[i];
				for (j=0; j<sb.Ribbons.Length; j++)
				{
					if( sb.Painter.RibbonRects[j].Contains(pbRibbonPt) )
					{
						Ribbon rib = sb.Ribbons[j];
						float fac = (pbRibbonPt.X-sb.Painter.RibbonRects[j].Left)/sb.Painter.RibbonRects[j].Width;
						return rib.Spine[Math.Min(rib.Spine.Length-1,(int)(fac*(float)rib.Spine.Length+0.5f))];
					}
				}
				for (j=0; j<sb.Bridges.Length; j++)
				{
					if( sb.Painter.BridgeRects[j].Contains(pbRibbonPt) )
					{
						OcclusionBridge bridge = sb.Bridges[j];
						float fac = (pbRibbonPt.X-sb.Painter.BridgeRects[j].Left)/sb.Painter.BridgeRects[j].Width;
						return new PointF( fac*bridge.Vf2.X+(1.0f-fac)*bridge.Vf1.X, fac*bridge.Vf2.Y+(1.0f-fac)*bridge.Vf1.Y );
					}
				}
			}
			//	Nothing found.
			return PointF.Empty;
		}
		public SurfacesDrawer Painter;
		public class SurfacesDrawer
		{
			public static Font InvalidFont = new Font("Arial", 20.0f);
			public static Brush InvalidBrush = new SolidBrush(Color.Red);
			#region Natural <--> View Transforms
			public Size PbSize
			{
				get
				{
					return pbSize;
				}
				set
				{
					pbSize = value;
					if(pbSize.Height>0.0f && pbSize.Width>0.0f)
						pbWhRatio=pbSize.Width/pbSize.Height;
					else
						pbWhRatio=Single.NaN;
					//	Reconfigure viewing rectangles.
					ViewerRect = viewerRect;
				}
			}
			private Size pbSize;
			private float pbWhRatio, natWhRatio;
			public RectangleF ViewerRect
			{
				get
				{
					return viewerRect;
				}
				set
				{
					viewerRect=value;

					if(Surfaces!=null)
					{
						this.CalcViewerRects();
						this.PassViewerRects();
					}
				}
			}
			private RectangleF viewerRect;
			private float l,r,t,b;
			/// <summary>
			/// Transforms a point from the viewerRect (painting area) to natural coordinates.
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float ViewerToNatX(float x)
			{
				return (x-viewerRect.Left) * PbSize.Width / (float)viewerRect.Width;
			}
			/// <summary>
			/// Transforms a point from the viewerRect (painting area) to natural coordinates.
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float ViewerToNatY(float y)
			{
				return (y-viewerRect.Top) * PbSize.Width / (float)viewerRect.Width;
			}
			public PointF ViewerToNat(PointF pt)
			{
				return new PointF(ViewerToNatX(pt.X), ViewerToNatY(pt.Y));
			}
			public RectangleF ViewerToNat(Rectangle rect)
			{
				l = ViewerToNatX((float)rect.Left);
				r = ViewerToNatX((float)rect.Right);
				t = ViewerToNatY((float)rect.Top);
				b = ViewerToNatY((float)rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			public RectangleF ViewerToNat(RectangleF rect)
			{
				l = ViewerToNatX(rect.Left);
				r = ViewerToNatX(rect.Right);
				t = ViewerToNatY(rect.Top);
				b = ViewerToNatY(rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			/// <summary>
			/// Transforms a point from natural coordinates to the viewerRect (painting area).
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float NatToViewerX(float x)
			{
				return x * (float)viewerRect.Width / PbSize.Width + viewerRect.Left;
			}
			/// <summary>
			/// Transforms a point from natural coordinates to the viewerRect (painting area).
			/// </summary>
			/// <param name="x"></param>
			/// <returns></returns>
			public float NatToViewerY(float y)
			{
				return y * (float)viewerRect.Width / PbSize.Width + viewerRect.Top;
			}
			public PointF NatToViewer(PointF pt)
			{
				return new PointF(NatToViewerX(pt.X), NatToViewerY(pt.Y));
			}	
			public RectangleF NatToViewerF(RectangleF rect)
			{
				l = NatToViewerX(rect.Left);
				r = NatToViewerX(rect.Right);
				t = NatToViewerY(rect.Top);
				b = NatToViewerY(rect.Bottom);
				return new RectangleF(l,t,r-l,b-t);
			}
			public Rectangle NatToViewer(RectangleF rect)
			{
				l = NatToViewerX(rect.Left);
				r = NatToViewerX(rect.Right);
				t = NatToViewerY(rect.Top);
				b = NatToViewerY(rect.Bottom);
				return new Rectangle((int)(l+0.5f),(int)(t+0.5f),(int)(r-l+0.5f),(int)(b-t+0.5f));
			}
			#endregion
			public SurfaceBoundList Surfaces;
			public RectangleF NaturalRect
			{
				get
				{
					if(SurfRects==null || SurfRects.Length==0)
						return RectangleF.Empty;
					else
						return new RectangleF(0.0f,0.0f,
							SurfRects[SurfRects.Length-1].Right, SurfRects[SurfRects.Length-1].Bottom);
				}
			}
			/// <summary>
			/// The natural rectangles for each surface.
			/// </summary>
			public RectangleF[] SurfRects, SurfRectsViewer;
			/// <summary>
			/// The natural rectangles for each empty space between surfaces.
			/// </summary>
			public RectangleF[] TweenRects, TweenRectsViewer;
			public float TweenWidth=10.0f;
			public Brush TweenBrush = new SolidBrush(Color.Firebrick);
			public void Paint(object sender, PaintEventArgs e)
			{
				if(Surfaces!=null && Surfaces.Count>0)
				{
					int i;
					//	First, fill all TweenRects with the TweenBrush
					for(i=0; i<TweenRects.Length; i++)
						e.Graphics.FillRectangle(TweenBrush, TweenRectsViewer[i]);
					//	Next, paint the surfaces
					for(i=0; i<Surfaces.Count; i++)
						Surfaces[i].Painter.Paint(sender, e);
					if(!Surfaces.IsValid)
						e.Graphics.DrawString("WARNING:  Invalid surfaces were specified", InvalidFont, InvalidBrush,
							this.viewerRect.X, this.viewerRect.Y);
					//	Next, call a ribbon to print the labels.
					Ribbon rib = AnyRibbon;
					if(rib!=null && TweenRects!=null && TweenRects.Length>0)
						rib.Painter.DrawText(sender, e, new PointF(
							Math.Min(SurfRectsViewer[0].Left, TweenRectsViewer[0].Left)-10.0f, SurfRectsViewer[0].Top ),
							SurfRectsViewer[0].Width/SurfRects[0].Width );
					else
						rib.Painter.DrawText(sender, e, new PointF(
							SurfRectsViewer[0].Left-10.0f, SurfRectsViewer[0].Top ),
							SurfRectsViewer[0].Width/SurfRects[0].Width);
				}
				else
				{
					e.Graphics.DrawString("The active object does not own any edges.",
						InvalidFont, InvalidBrush, 5.0f, 5.0f, StringFormat.GenericTypographic);
				}
			}
			public Ribbon AnyRibbon
			{
				get
				{
					if(Surfaces!=null && Surfaces.Count>0)
					{
						for(int i=0; i<Surfaces.Count; i++)
						{
							if(Surfaces[i].Ribbons!=null && Surfaces[i].Ribbons.Length>0)
								return Surfaces[i].Ribbons[0];
						}
					}
					return null;
				}
			}
			public void CalcNaturalRects()
			{
				if(Surfaces!=null && Surfaces.Count>0)
				{
					//	Create a surface rect for each Surfaces item.
					SurfRects = new RectangleF[Surfaces.Count];
					TweenRects = new RectangleF[Surfaces.Count-1];
					SurfRectsViewer = new RectangleF[Surfaces.Count];
					TweenRectsViewer = new RectangleF[Surfaces.Count-1];
					int i;
					float left=5.0f;
					//	For each surface rect (Surfaces item)
					for(i=0; i<SurfRects.Length; i++)
					{
						//	Do Tweener that comes before this surface rectangle.
						if(i>=1)
						{
							TweenRects[i-1] = new RectangleF(
								SurfRects[i-1].Right, SurfRects[i-1].Top,
								TweenWidth, SurfRects[i-1].Height);
						}
						//	Create the surface rectangle.
						SurfRects[i] = new RectangleF(left, 0.0f,
							Surfaces[i].Painter.NaturalWidth, Surfaces[i].Painter.NaturalHeight);
						//	Setting this property item initiates a cascade of changes in the lower objects
						//	whereby all the natural rects are reset.
						Surfaces[i].Painter.NaturalRect = SurfRects[i];
						//	Increment the left point for next run.
						left=SurfRects[i].Right+TweenWidth;
					}
				}
				//	Set natural ratio
				RectangleF nat = NaturalRect;
				natWhRatio=nat.Width/nat.Height;
			}
			public void CalcViewerRects()
			{
				if(Surfaces!=null && Surfaces.Count>0)
				{
					int i;
					//	For each surface rect (Surfaces item)
					for(i=0; i<SurfRects.Length; i++)
					{
						//	Do Tweener that comes before this surface rectangle.
						if(i>=1)
							TweenRectsViewer[i-1] = this.NatToViewer(TweenRects[i-1]);
						//	Create the view rects
						SurfRectsViewer[i] = this.NatToViewer(SurfRects[i]);
					}
				}
			}
			public void PassViewerRects()
			{
				if(Surfaces!=null && Surfaces.Count>0)
				{
					int i;
					//	For each surface rect (Surfaces item)
					for(i=0; i<SurfRects.Length; i++)
					{
						//	Setting this property item initiates a cascade of changes in the lower objects
						//	whereby all the viewer rects are reset.
						Surfaces[i].Painter.ViewerRect=SurfRectsViewer[i];
					}
				}
			}
			/// <summary>
			/// Initializes a SurfacesDrawer object.  Ribbons must already be constructed.
			/// </summary>
			/// <param name="PbSize">The size of the picture-box control.</param>
			/// <param name="list">The list of surface bounds upon which this object is to be based.</param>
			public SurfacesDrawer(Size PbSize, SurfaceBoundList Surfaces)
			{
				this.Surfaces=Surfaces;
				this.pbSize = PbSize;
				CalcNaturalRects();
				//	Default viewer rects (initiate viewer cascade).
				ViewerRect = new RectangleF(new PointF(0.0f,0.0f), PbSize);
			}
		}
		#endregion
		#region CollectionBase
		public SurfaceBound this[int index]
		{
			get  
			{
				return( (SurfaceBound) List[index] );
			}
			set  
			{
				List[index] = value;
			}
		}
		public int Add( SurfaceBound sb )
		{
			return List.Add( sb );
		}

		public int IndexOf( SurfaceBound sb )
		{
			return( List.IndexOf( sb ) );
		}

		public void Insert( int index, SurfaceBound sb )  
		{
			List.Insert( index, sb );
		}

		public void Remove( SurfaceBound sb )  
		{
			List.Remove( sb );
		}
		public new void RemoveAt(int index)
		{
			List.RemoveAt( index );
		}
		public bool Contains( SurfaceBound sb )  
		{
			return List.Contains( sb );
		}
		#endregion
	}
}